CONTRACT  REPORT  NO.  334 

NON  LINEAR  DYNAMIC  ANALYSIS  OF  FLAT 


\D 


LAMINATED  PLATES  BY  THE  FINITE-ELEMENT 


METHOD 


Prepared  by 

University  of  Illinois 
Aeronautical  and  Astronautical 
Engineering  Department 
Urbana,  Illinois 


STINFO  BRANCH 
BRL,  APG,  HD.  21005 


March  1977 


> 


Approved  for  public  release;  distribution  unlimited. 


USA  ARMAMENT  RESEARCH  AND  DEVELOPMENT  COMMAND 
USA  BALLISTIC  RESEARCH  LABORATORY 

ABERDEEN  PROVING  GROUND,  MARYLAND 


Destroy  this  report  when  it  is  no  longer  needed. 
Do  not  return  it  to  the  originator. 


Secondary  distribution  of  this  report  by  originating 
or  sponsoring  activity  is  prohibited. 

Additional  copies  of  this  report  may  be  obtained 
from  the  National  Technical  Information  Service, 

U.S.  Department  of  Commerce,  Springfield,  Virginia 
22151. 


The  findings  in  this  report  are  not  to  be  construed  as 
an  official  Department  of  the  Army  position,  unless 
so  designated  by  other  authorized  documents. 


SECURITY  CLASSIFICATION  OF  This  PAGE  (When  Data  Entered) 


REPORT  DOCUMENTATION  PAGE 


1.  REPORT  NUMBER 


lontract  Report  Number  334' 


4.  TITLE  (and  Subtitle) 

Non  Linear  Dynamic  Analysis  of  Flat  Laminated 
Plates  by  the  Finite- Element  Method 


9.  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

Aeronautical  and  Astronaut ical  Engr.  Dept 
University  of  Illinois 
Urbana,  Illinois 


11.  CONTROLLING  OFFICE  NAME  AND  ADDRESS 

USA  Ballistic  Research  Laboratory 
Aberdeen  Proving  Ground,  Maryland  21005 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


3.  RECIPIENT'S  CATALOG  NUMBER 


TYPE  OF  REPORT  & PERIOD  COVERED 

BRL  Contractor  Report 
1 Jan  75  to  30  Jun  76 


PERFORMING  ORG.  REPORT  NUMBER 

None 


CONTRACT  OR  GRANT  NUMBERf*} 


DAAD05-73-C-0197 


10.  PROGRAM  ELEMENT.  PROJECT,  TASK 
AREA  ft  WORK  UNIT  NUMBERS 

1T161702A33H,  002AJ 
61702A 


12.  REPORT  DATE 

MARCH  1977 


13.  NUMBER  OF  PAGES 


4.  MONITORING  AGENCY  NAME  ft  ADDRESSflf  different  from  Controlling  Office)  15.  SECURITY  CLASS,  (of  thlm  report) 

US  Army  Materiel  S Readiness  Command 

5001  Eisenhower  Avenue  Unclassified 

Alexandria,  VA  22333 


15«.  DECLASSIFICATION/DOWNGRADING 
SCHEDULE 


16.  DISTRIBUTION  STATEMENT  (of  thtm  Report) 


pproved  for  public  release;  distribution  unlimited. 


17.  DISTRIBUTION  STATEMENT  (of  the  abstract  entered  in  Block  20,  It  different  from  Report) 


16.  SUPPLEMENTARY  NOTES 


PROPERTY  OP  U.S.  ARMY 
S.TINF0  BRANCH 
BRL,_  APG,  MQ..  21005 


19.  KEY  WORDS  (Continue  on  rover  me  eidm  if  noce  emery  end  identify  by  block  number) 

Anisotropic,  Analysis,  Finite-Element,  Three  Dimensional  Deformations, 
Dynamic  Response 


2CL  ABSTRACT  (Contfmrm  mm  wo 


l Identity  by  block  number) 


A finite-element  structural  model  has  been  developed  for  the  dynamic  analysis 
of  laminated,  thick  plates.  The  model  uses  constant  thickness  quadrilateral 
elements  to  represent  the  shape  of  the  plate  and  these  elements  are  stacked 
in  the  thickness  direction  to  represent  the  desired  material  layers.  The 
analysis  allows  for  orthotropic  material  properties  of  each  layer  as  well  as 
for  elastic-plastic  material  response.  Nonlinear  strain  displacement  relations 
are  used  in  the  formulation  to  represent  large,  transverse  plate  deflections. 


FORM 
I JAM  73 


EDITION  OF  V NOV  65  IS  OBSOLETE 


UNCLASSIFIED 

. SECURITY  CLASSIFICATION  OF  THIS  PAGE  (When  Date  Entered) 


SECURITY  CLASSIFICATION  OF  THIS  PAOE{UTi«  Datm  Bntmd) 


The  finite-element  model  is  used  to  prepare  a computer  program  for  the 
numerical  calculations.  ' Two  versions  of  the  program  have  been  prepared, 
which  correspond  to  two  different  time  integration  numerical  methods.  These 
methods  include  finite-difference  and  predictor- corrector  techniques.  The 
computer  programs  are  designed  for  time  and  space  dependent  pressure  loads 
to  be  applied  to  one  surface  of  the  plate.  However,  the  programs  could  be 
used  for  other  loading  conditions,  by  changing  one  subroutine. 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  PAGEfWIwi  Data  E«f«r®0 


TABLE  OF  CONTENTS 


Page 

I.  INTRODUCTION 5 

II.  THEORETICAL  DEVELOPMENT . 6 

Finite-Element  Model  6 

Dynamic  Equations  7 

III.  NUMERICAL  ANALYSIS  11 

Finite-Difference  Method  11 

Predictor-Corrector  Method  .....  13 

Elastic-Plastic  Analysis  13 

Computer  Program  Inputs  15 

Solution  of  Laminated  Plate  16 

IV.  CONCLUSIONS  .....  17 

FIGURES . . 18 

APPENDIX  A DERIVATION  OF  DIFFERENCE  EQUATIONS  ....  27 

APPENDIX  B COMPUTER  PROGRAM  INPUT  CARD  . 31 

APPENDIX  C COMPUTER  PROGRAM  LISTING  DESCRIPTION 

FOR  FINITE  DIFFERENCE  PROGRAM  39 

DISTRIBUTION  LIST 81 


3 


I INTRODUCTION 


The  purpose  of  this  investigation  has  been  to  develop  a finite- 
element  model  and  a computer  program  for  the  dynamic  analysis  of  flat, 
laminated  plate  structures.  The  finite-element  model  uses  the  quadrila- 
teral to  define  the  shape  of  the  element  in  the  plane  of  the  plate  and 
the  thickness  direction  is  represented  by  arranging  a number  of  these 
elements  to  describe  the  necessary  number  of  material  layers.  Each 
material  layer  can  be  assigned  different  material  properties  and  the 
computer  program  is  designed  to  allow  for  maximum  of  twelve  layers, 
however,  this  can  easily  be  increased  or  decreased.  In  the  detail  de- 
velopment of  the  model  each  quadrilateral  element  is  further  subdivided 
into  four  triangular  elements! >2. 

The  model  in  this  analysis  is  nonlinear  since  it  allows  for  material 
yield  effects  and  for  large  plate  deflections . The  large  deflection  ef- 
fects are  introduced  by  assuming  that  the  transverse  displacement  is 
large  compared  to  the  two  displacements  in  the  plane  of  the  plate.  This 
leads  to  second  order  terms  in  the  strain-displacement  relations.  The 
yield  effects  are  introduced  because  it  is  intended  that  the  application 
of  this  model  will  be  in  the  situation  where  the  loads  are  large  enough 
to  produce  stresses  beyond  the  elastic  limit.  Consequently,  the  present 
analysis  allows  for  elastic-plastic  material  properties  which  are  intro- 
duced into  the  model  by  checking  the  yield  for  each  element  and  when  the 
yield  is  exceeded,  the  state  of  stress  in  the  element  is  adjusted  by  using 
the  plasticity  flow  rule.  This  is  done  by  checking  for  yield  at  each  time 
interval  used  in  the  numerical  integration  of  the  dynamic  equations. 

The  dynamic  equations  for  the  plate  are  obtained  by  lumping  the 
mass  of  the  plate  into  the  nodal  points  of  the  finite-element  model. 

This  leads  to  a set  of  concentrated  masses  distributed  in  the  plane  of 
the  plate  and  in  the  thickness  direction.  The  solution  to  these  equations 
is  obtained  numerically  in  the  computer  program.  In  the  preliminary  ver- 
sion of  the  program  three  different  integration  techniques  are  investigated. 
These  include  an  iterative  approach,  a finite-difference  method,  and 
a predictor-corrector  method.  Based  on  some  results  obtained  for  simple 
dynamic  problems,  it  was  found  that  the  fastest  methods  were  the  finite- 
difference  and  the  predictor-corrector  with  the  iterative  approach  being 


^Zienkiewicz,  0.  C.,  "The  Finite  Element  Method  in  Engineering 
Science,"  McGraw-Hill  Publishing  Company,  London,  1971. 
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Przemieniecki,  J.  S.,  "Theory  of  Matrix  Structural  Analysis," 
McGraw-Hill  Book  Company,  New  York,  1968. 


the  slowest.  This  of  course  does  not  imply  that  the  same  results  would 
be  obtained  for  all  dynamic  problems,  however,  in  the  present  investi- 
gation it  was  necessary  to  limit  the  number  of  versions  of  the  computer 
program  and,  therefore,  it  was  decided  to  prepare  two  final  versions 
using  the  finite-difference  and  the  predictor-corrector  methods  of 
integration. 

The  analysis  which  has  been  developed  is  quite  general  and  corres- 
ponding computer  program  could  be  used  for  a wide  class  of  problems 
and  loading  conditions.  Various  loading  conditions  can  be  generated 
by  supplying  user  subroutine  to  define  the  load  in  space  and  time. 
However,  the  present  versions  of  the  program  is  setup  for  particular 
load  which  involves  a distributed  pressure  bn  one  surface  of  the  plate. 
This  pressure  can  be  specified  as  a function  of  time  and  of  the  inplane 
coordinates. 


II  THEORETICAL  DEVELOPMENT 

Finite-Element  Model 

The  present  analysis  is  intended  to  handle  thick,  plate-like 
structures,  which  are  composed  of  different  material  layers  in  the 
thickness  directions.  The  plane  of  the  plate  is  parallel  to  the  x-y 
coordinates  and  the  thickness  is  represented  by  the  z direction  as  shown 
in  Figure  1.  The  analysis  is  developed  in  terms  of  these  Cartesian  co- 
ordinates. The  shape  of  the  finite  elements  is  defined  by  a general 
quadrilateral  in  the  x-y  plane  and  each  element  has  a constant  thickness 
in  the  z direction  with  the  limitation  that  no  element  will  contain  more 
than  one  material.  Each  layer  of  the  material  can  therefore  be  repre- 
sented by  one  or  more  elements  in  the  thickness  direction. 

A typical  element  is  shown  in  Figure  2.  Each  quadrilateral  element 
is  subdivided  into  four  triangular  elements  with  the  node  in  the  center 
being  a temporary  that  will  later  be  eliminated  by  static  condensation. 

A general  triangular  element  is  shown  in  Figure  3.  The  nodal  numbering 
system  for  each  quadrilateral  element  is  also  given.  There  are  three 
degrees  of  freedom  at  each  node. 

The  first  step  in  the  finite  element  analysis  is  to  assume  a suitable 
displacement  function  over  each  triangular  element.  The  displacement  func 
tions  chosen  for  this  analysis  are 

0=0^+  a2x  + a3y  + z (Bj  + B2x  + B3y) 

v = a4  + agx  + a6y  + z (34+  35x  + B6y) 

w = a?  + agx  + agy  + z (3?+  3gx  + 3gy)  (1) 

where  a.,  i = 1,  9,  and  6.,  i = 1,  9,  are  unknown  coefficients.  This 
displacement  function  is  linear  in  planes  parallel  to  the  x-y  plane  and 
varies  linearly  with  z through  the  thickness  of  each  element.  Writing 
Equations  (1)  at  each  of  the  six  nodes  of  a triangular  element  results 
in  eighteen  equations  which  can  be  written  in  matrix  form  as 


IV 


{6}  = [C]  {a} 


(2) 


where  {6}  denotes  the  nodal  displacements,  [C]  is  a known  constant  mat- 
rix depending  on  the  local  nodal  coordinates,  and  {a}  is  the  vector  of 
unknown  coefficient  a.  and  8 defined  in  Equations  (1).  Solving  for  {a} 
results  in 


{a}  = [C]"1  {6} 


(3) 


It  may  be  noted  at  this  time  that  as  a consequence  of  choosing  the  dis- 
placement variations  as  given  by  Equations  (1),  the  inverse  of  the  [C] 
matrix  in  equation  (3)  can  be  performed  analytically  thereby  leading  to 
an  appreciable  saving  in  numerical  work.  The  present  analysis  is  de- 
signed to  account  for  nonlinear  effects  arising  from  large  deflections 
of  the  plate.  In  order  to  account  for  these  deflections  analytically, 
it  is  assumed  that  the  deflections  and  rotations  out  of  the  plane  of 
the  plate  are  larger  compared  to  those  in  the  plane  of  the  plate.  The 
resulting  non-linear  strain-displacement  relations  are  therefore  given  by 
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Equations  (4)  may  be  written  in  term  of  linear  and  nonlinear  terms.  Using 
matrix  notation 


where 


{e}  = {e  } + {e  } 
0 L 
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(6) 


(7) 


Dynamic  Equations 

The  next  step  in  the  finite  element  model  derivation  is  the  calcu- 
lation of  the  virtual  change  in  the  internal  work  of  the  structure  due  to 
virtual  changes  in  the  nodal  displacements.  The  virtual  change  in  the  in- 
ternal work  due  to  the  stresses  is  given  by  an  integral  over  the  volume  of 
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a particular  element.  This  can  be  represented  by 


dW, 


T • 1 

d{e}  {a}  dV 
■ V 


(8) 


T 

where  dWT  represents  the  virtual  change  in  the  internal  work,  d{e}  is  the 
transpose  of  the  matrix  representing  the  virgual  changes  in  the  strains, 
and  {a}  is  the  matrix  of  the  stresses.  The  quantity  V represents  the  volume 
of  the  element.  The  first  step  in  obtaining  the  integral  in  Equation  (8) 
is  to  derive  the  expression  for  the  virtual  changes  in  the  strains.  Using 
Equation  (5)  written  in  two  parts 


d{e}  = d{£  } + d{er } 

O L 


(9) 


Consider  first  the  linear  strain  component  of  virtual  strain  d{eQ}.  By 
using  the  definition  of  {e  } given  by  Equation  (6)  and  the  displacement 
function  in  Equations  (1)  it  is  possible  to  write 


{eQ>  = [Q]  ia}  (10) 

where  [Q]  is  a matrix  whose  elements  are  functions  of  x,  y and.z.  By 

combining  Equations  (3)  and  (10)  it  follows 

(S0)  . [1?]  «>'  (11) 

where  the  [B  1 matrix  is  defined  as 

o 

[B0]  = [Q]  [C]'1  (12) 


By  introducing  the  virtual  change  of  nodal  displacements  it  follows  that 


d{e  } - [B  ] d{6)  (13) 

o o 

where  d{6}  are  changes  in  the  nodal  displacements. 


Consider  now  the  virtual  change  in  the  nonlinear  component  of  strain 

d{e  }.  Consider  Equation  (7)  which  can  be  rewritten  in  a different  form 
L 

{eL}  = \ [A]  {9}  (14) 

where  the  new  matrices  in  Equation  (14)  are  defined  as  follows 
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It  may  be  noted  that  using  Equations  (15)  the  strain  {e.}  has  only  four 
elements  while  the  original  definition  of  in  Equation  (5),  has  six 

elements.  However,  it  may  be  noted  from  the  original  definition  that  the 
last  two  elements  {e,}  were  identically  zero  and  consequently  (e. } in 
Equation  (13)  represents  only  the  non  zero  part  of  the  original  strains. 
The  virtual  change  in  the  strain  {eL)  can  be  written  as 

dUL>  = i d[A]  {0}  ♦ j [A]  d{9)  (16) 

However,  by  using  Equations  (15)  it  can  be  shown  that 


d[A]  {9}  = [A]  d{6) 

therefore 


(17) 


d(eL)  = [A]  d{0)  (18) 

By  using  the  definition  of  d{0)  it  is  possible  to  write 

d{0)  = [X]  d{6}  (19) 

It  may  be  noted  at  this  stage  that  the  matrix  [A]  contains  the  coordi- 
nates x,  y,  z and  some  of  the  unknown  coefficient  a.  and  0.  defined  in 
Equations  (1) . During  the  programming  of  this  analysis  for  numerical  cal- 
culations, it  was  found  convenient  to  write  Equation  (19)  in  slightly  dif- 
ferent form 


d(eL)  = [Z]  [C]  d{6)  (20) 

where  obviously 

[Z]  [C]  = [A]  [X]  (21) 

The  matrix  [C]  in  Equation  (21)  is  actually  related  to  the  matrix  [C]~ 
defined  in  Equation  (3)  and  [Z]  contains  the  variables  x,  y,  z and  the 
cl  and  0^  coefficients. 

Returning  to  Equation  (8) , the  internal  work  term  can  be  written  as 
dWT  = d{6}T  f [B  ]T  {a}  dV 

I JV  ° 

+ d(6}T  [C]T  [ [Z]T  {a}  dV  (22) 

JV 

The  stress  matrix  {a}  is  defined  as  follows 

(a)  * [ on,  o22,  a33,  a12,  a23,  a3J]  (23) 


However,  it  may  be  noted  that  in  the  second  integral  in  Equation  (22)  only 
the  first  four  stresses  from  Equation  (23)  are  needed.  By  combining  the 
two  integrals  in  Equation  (22)  it  is  possible  to  write 
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(24) 


6WX  = d{<S}T 


[B] 


T {cr}  dV 


V 


where  {a}  is  defined  according  to  Equation  (23) . 


In  addition  to  the  work  done  by  internal  forces  in  each  finite  ele- 
ment, Equation  (24),  there  is  work  done  by  the  forces  at  the  eight  corners. 
By  denoting  these  forces  by  the  matrix  {f},  the  external  virtual  work 
done  is  obtained  in  the  form 

dW_  = d(6}T  {f } (25) 

E 


Since  the  external  forces  on  each  element,  are  in  equilibrium  with  the  stresses 
within  the  element  it  follows  from  the  principle  of  virtual  work  that 


6Wj  + 6WE  = 0 


and  therefore 


(f)  = - 


] [b]t  {a}  dV 

Jv 


(26) 


(27) 


The  next  step  in  the  analysis  is  to  obtain  the  matrix  equilibrium 
equation  for  the  total  structure.  Before  this  can  be  done,  it  is  necessary 
to  introduce  the  inertia  effects.  In  the  present  analysis  the  mass  of  the 
structure  is  lumped  at  the  nodes.  Each  side  of  the  triangular  element 
shown  in  Figure  3 is  bisected  and  these  center  points  are  joined  to  the. 
center  node  which  defines  the  apex  of  the  four  triangular  elements  previously 
defined  in  Figure  2.  By  performing  this  step  the  original  quadrilateral  has 
been  divided  into  four  smaller  quadrilaterals  and  each  of  these  contains  two 
of  the  eight  nodes  of  the  element.  The  mass  of  each  smaller  quadrilateral 
portion  of  the  element  is  obtained  by  multiplying  the  area  of  the  quadri- 
lateral by  the  thickness  of  the  element  in  the  z direction  and  by  the  ma- 
terial density  and  this  mass  is  distributed  equally  to  the  two  nodes. 
Following  this  procedure,  the  mass  matrix  for  the  total  structure  is  ob- 
tained by  summing  masses  at  each  structure  node  from  the  adjacent  elements. 
The  resulting  matrix  is  diagonal. 

The  next  step  of  the  analysis  is  the  summation  of  forces  at  each  node 
of  the  structure.  The  forces  acting  on  any  node  from  the  adjacent  elements 
are  obtained  from  Equation  (27) . The  sum  of  these  forces  can  be  represented 
for  the  total  structure  by  a matrix  which  will  be  denoted  by  (Fj);  The 
additional  forces  which  affect  the  system  are  the  external  forces.  The  pre- 
sent model  will  allow  for  these  to  be  surface  pressure  forces,  which  vary 
with  space  and  time.  Using  the  principle  of  virtual  work  these  pressure 
loads  are  transformed  to  an  equivalent  set  of  concentrated  nodal  forces. 

These  forces  will  be  denoted  by  a matrix  {F£}.  Using  the  two  forces  the 
matrix  equilibrium  equation  can  be  written  in  the  form 


[M]  {&}  = (Fj)  + {Fe'}  (28) 

where  [M]  is  the  mass  matrix,  and  (A)  represents  the  global  displacement 
matrix  with  the  double  dot  denoting  the  second  time  derivative.  In  the 
present  investigation  the  solution  to  Equation  (28)  is  obtained  using 
numerical  integration  methods  which  will  be  described  in  the  next  section 


<• 
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of  this  report. 


I I I NUMERICAL  ANALYS I S 


The  solution  of  Equation  (28)  is  obtained  by  using  two  different 
numerical  integration  techniques.  These  two  different  methods  include  a 
finite-difference  approach  and  a predictor-corrector  method.  Both  of  these 
methods  were  used  to  prepare  two  different  versions  of  the  computer  pro- 
gram. A brief  description  of  both  of  these  methods  is  given  below. 

Finite  Difference  Method 


The  numerical  solution  of  Equation  (28).  by  the  finite  difference 
method  involves  the  replacement  of  the  time  derivatives  by  finite  difference 
equivalents  3,4,5,  in  this  approach  the  time  history  is  .divided  into  dis- 
creet intervals  whose  length  will  be  denoted  by  h.  For  convenience  any 
general  element  of  the  displacement  vector  {A}  at  a given  time  will  be 
represented  by  the  quantity  x with  a subscript  defining  the  time  interval . 
From  the  general  theory  of  kinematics,  the  velocity  and  displacement  re- 
lations are  written  for  nth  and  (n+1)  th  time  intervals  as  follows 


* * h ...  ..  . 

X = X + — (x  - + X ) 

n n-1  2 v n-1  n' 


x = x . + hx  . + (-iy  - 8)h2  x 1 + Bh2x 

n n-1  n-1  v2  n-1  n 


• * h ..  h .. 

x = x + •=■  x + ■?  x 
n+1  n 2 n 2 n+1 


x ,=x  + h x +(4-  -8)h2x  + Bh2x 

~ " n 2 n 


n+1  n 


h+1 


(29) 

(30) 

(31) 

(32) 


* 

where  x,  x,  and  x represent  acceleration,  velocity  and  displacements  at 
the  time  intervals  denoted  by  the  subscript. 


Equations  (29)  and  (31)  mean  that  the  velocity  at  the  end  of  the  in- 
terval is  equal  to  the  sum  of  the  velocity  at  the  beginning  of  the  interval 
and  the  product  of  the  length  of  the  interval,  h,  and  the  average  of  the 
accelerations  at  the  beginning  and  end  of  the  interval . Equations  (30) 
and  (32)  are  obtained  by  integrating  Equations  (29)  and  (31)  and  introducing 


3 

Newmark,  Nathan,  M. , "A  Method  of  Computation  for  Structural  Dynamics," 
Journal  of  the  Engineering  Mechanics  Division,  ASCE,  July,  1959. 

4 

Chan,  S.  P.,  Cox,  H.  L.,  and  Benfield,  W.  A.,  "Transient  Analysis 
of  Forced  Vibrations  of  Complex  Structural-Mechanical  Systems,"  Journal 
of  the  Royal  Aeronautical  Society,  July,  1962. 

5 

Wu,  R. , and  Witmer,  E.  A.,  "Nonlinear  Transient  Responses  of 
Structures  by  the  Spatial  Finite  Element  Method,"  AIAA  J. , August,  1973. 
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the  acceleration  parameter,  B,  to  express,  the  acceleration  at-  the  be- 
ginning and  end  of  the  interval.  For  example,  3 = j means  that  the 
acceleration  during  the  interval  is  constant  and  is  equal  to  the  mean  . 
of  the  accelerations  at  the  beginning  and  end  of  the  time  interval. 

The  3 parameter  is  known  as  the  generalized  acceleration 
parameter  and  its  value  is  chosen  in  the  numerical  calculations  to 
insure  convergence  and  stability  of  the  numerical  results.  According 
to  previously  published  results'* it  has  been  found  that  in  order  to  insure 
stability  and  convergence  the  value  of  this  parameter  should  be  kept  in 
the  range  of  0 < B < 1/4.  The  different  values  of  this  parameter  are 
suitable  for  different  types  of  dynamic  problems,  however,  when  B = 1/4  » 

the  stability  limit  on  the  integration  time  interval  is  infinite.  Con- 
sequently, this  value  of  3 should  be  suitable  in  most  dynamic  problems.’ 

« 

At  t = (n+1)  h,  nh,  and  (n-l)h,  respectively,  the  equation  of  motion, 

Equation  (28),  becomes 

[M]  tonil  - (Fj)  - <Fe) 

[«]  <*>„  - tFI>  n - {FE>  n t33) 

f]  - tFli  n-1  * (FE>  n-1 

By  combining  Equations  (29)  to  (33)  it  can  be  shown,  see  Appendix  A,  that 
this  leads  to  the  following  difference  equation  for  the  displacement  at 
the  n+1  time  interval 

[M]U)ntl  = 2[M]{A>n  - [M]{A)n  l 

* Bh2‘{FI>  n+1  * (1/6  - 2){Fl!  n * {Fli  n-13  <34) 

♦ 61,2 ((Fe3  n+1  * 2)(Fe)  n * (Fe)  n.p 

Equation  (34)  is  used  in  the  computer  program  to  calculate  the  displace- 
ment at  the  time  interval  n+1  from  the  response  at  the  two  previous  time 
intervals  n and  n-1.  In  applying  Equation  (34)  it  is  assumed  that  the  in- 
ternal force  (Ft)  - is  dependant  on  the  stresses  calculated  at  the  time 
interval  n.  The’pRysical  interpretation  is  that  the  analysis  treats  the 
internal  force  as  a step  function  rather  than  as  a continuous  force.  This  , 

can  be  illustrated  by  a curve  shown  in  Figure  4 , which  represents  a force 
constant  over  each  time  interval.  Obviously  this  is  an  approximation, 
however,  as  the  time  interval  decreases  this  curve  will  approach  a con- 
tinously  varying  function. 

It  may  be  noted  that  Equation  (34)  can  not  be  applied  to  the  first 
time  interval  since  {A}  , does  not  exist.  Consequently,  in  the  beginning 

of  this  solution  procedQre  the  displacement  at  the  end  of  the  first  inter- 
val must  be  formulated  in  terms  of  the  initial  velocity  and  initial 
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displacement  since  these  are  the  only  known  quantities.  It  is  possible  to 
use  the  following  starting  procedure 

[Ml  {Alj  + 0h2  (-{FjJj)  = [M]  {A}q 

- (j-  6)h2  (-{Fj>  q)+  [M]h  {A} 

♦ Bh2  (FE>  j + (I  - 6) [ I ] h2  {Fe>  o . (35) 

where  [I]  is  the  unit  matrix.  Equation  (35)  is  derived  in  a manner  similar 
to  Equation  (34) . 

Thus,  the  general  procedure  is  to  start  with  Equation  (35)  in  order 
to  obtain  the  displacements  at  the  end  of  the  first  time  interval  and 
subsequently  use  Equation  (34)  to  obtain  the  displacements  at  later  times. 
(Note  that  when  8=0,  this  method  reduces  to  the  central  finite  difference 
method. ) 

Predictor-Corrector  Method 

In  general,  predictor-corrector  methods  involve  using  a truncated  for- 
mula to  'predict'  the  value  of  the  unknown  and  then  applying  a more  accurate 
'corrector'  formula  to  provide  successive  improvements. 

The  predictor-corrector  subroutine  used  in  this  analysis  is  named  DHPCG 
(Hamming's  Predictor-Corrector  Method)  and  is  from  the  IBM-  Scientific  Sub- 
routine Package  which  also  gives  a detailed  explanation  of  the  procedure  . 

In  brief,  Hamming's  Predictor-Corrector  method  gives  an  approximate  numeri- 
cal solution  to  a first  order  linear  ordinary  differential  equation  with 
given  initial  conditions.  It  is  a stable  fourth  order  procedure  in  which 
the  user  may  vary  the  step-size.  DHPCG  also  estimates  the  local  truncation 
error. 

Elastic-Plastic  Analysis 

The  present  analysis  allows  for  elastic-plastic  response  of  the  struc- 
ture. Because  of  practical  considerations  each  finite  element  is  assumed 
to  be  either  elastic  or  plastic.  However,  the  stress  strain  relation  which 
can  be  written  in  a matrix  form  as 

(a)  = [D]  {e)  (36) 


^Ralston,  A.,  and  Wilf,  H.  S. , "Mathematical  Methods  for  Digital 
Computers,"  Wiley,  New  York,  1960,  p.  95-109; 


predicts  that  the  stresses,  will  vary  within  each  element  since  {e},  as 
given  by  Equation  (5)  is'a  function  of  the  coordinates.  In  order  to  have 
one  representation  yield  criterion  for  each  element,  a numerical  proce- 
dure was  adopted,  which  is  used  in  Equation  (36),  which  averages  the. 
stresses  over  each  element.  . These  average  stresses  were  used  in  evaluating 
the  internal  nodal  forces  in  Equation  (27). 

The  Mises-Hencky  yield  .criterion5  is  used  to  determine  if  plastic 
flow  has  occured  in  any  given  element.  The  Stresses  and  strains,  which 
here  will  be  denoted  by  indicial  notation  at  time  t^  and  the  displace- 
ments at  time  t are  known. . From  this  information  the  stress  increment 
at  time  t . isncalculated  from  the  elastic  constitutive  relations  , 


that  is , 


(Aaij)Tn+l  = 1+v  I (Aeij)n+l  + l-2v 


kk'n+1  6ij 


(37). 


The  total  stress  at  t .is 

n+l 


(a,.)  = (a.,)  + (A  a.,) 


n+l 


ir 


n 


il' 


n+l 


where  the  superscript  ' T ' denotes  a trial  stress  state  and  (a4j)n+1  is 

total  trial  stress  at  time  t .,  (0..)  is  the  total  stress  atJtime  tn, 

and  (Aa. .)T  , is  the  trial  stress  increment  at  time  t . . 
v iyn+1  n+l 

The  Mises-Hencky  yield  function  is  given  by  • 


(38) 

the 


$ = 


S.  . 
il 


a 2 

y 


(39) 


where  S. . are  the  deviatoric  components  of  stress,  that  is, 
il 

S.  . - a.  . - y a , 6. 

3 kk  ij 


(40) 


and  a is  the  known  uniaxial  yield  stress  of  the  material. 

y 

Substituting  Equation  (37)  into  Equation  (38)  gives 

$T  = (S..)T  . (S..)T  . a 2 (41) 

n+l  v ij  n+l  il  n+l  3 y 

If  < 0 then  the  trial  stress  state  is  in  the  elastic  region  and  no 

plastic  flow  has  occured.  In  this  case  the  total  stress  is  given  by  Equa- 
tion (38)  or 


CO  ♦(*»,.) 


n+l 


il 


n 


i r 


n+l 


(42) 


If  $T  , ^ 0 then  plastic  flow  has  occured  and  the  stress  increments  is 
not  totally  elastic  as  was  assumed.  The  stress  state  must  lie  on  the  yield 
surface  as  specified  by'  the  theory  of  perfect  plasticity.  To  calculate  the 
new  stress  state,  the  strain  increment  is  broken  down  into  elastic  and 
plastic  components,  that  is. 


‘v  J I 

•J 


14 


(43) 


(Ae. .) 

±y 


n+1 


t'V 


n+1 


(Ae . . ) 

i y 


P 

n+1 


where  the  superscripts  ' e'  and  ' p ' denote  the  elastic  and  plastic  components, 
respectively.  From  the  incompressibility  condition  of  plasticity  and  by 
the  flow  rule,  the  plastic  strain  increment  is  given  by 


(Ae. .) 

iy 


n+1 


(S.  .) 
ir 


n+1 


where  X is  a real  nonnegative  scalar  quantity. 


(44) 


Combining  Equation  (44)  with  Equation  (37)  the  stress  increment  is 

.t  :] 


(Act.  .) 

ij 


n+1 


E rA  V 

1+v  1 ij  l-2v  “wkk  wij 


Ae,.,.  <5, 


(s.  .r  , XJ 

ij  n+1 


(45) 


and  the  actual  stress  at  time  t 


n+1 


(CT.  .) 

ij 


n+1 


= (CT..) 

ij 


n+1 


- (S.  .) 
ij 


X* 


(46) 


n+1 


where 

X*  = XE/l+v 

The  quantity  X*  may  be  determined  from  the  fact  that  (ct..)  1 in  Equation 

(46)  must  satisfy  the  yield  criterion  that  is  $ = 0.  Substituting  Equa- 
tion (46)  into  Equation  (39)  and  solving  for  X*  gives 


X*  - 


c 

B+/  B*  - AC 


(47) 


where  for  convenience  the  following  parameters  were  introduced 


T T 

a = cs,  ,y  (s,  .y 


n+1 


ij 


n+1 


T T 

8 = (°iP  ■ 

C * - T -J°y2‘  Vl 


(48) 


Computer  Program  Inputs 


The  computer  programs  for  the  two  different  versions  of  integration 
have  been  designed  to  have  the  same  type  of  input  cards.  The  description 
of  the  input  cards  is  given  in  Appendix  B.  The  two  versions  of  the  com- 
puter program  are  quite  similar  as  seen  from  the  flow  charts  shown  in 
Figures  5 and  6,  which  are  for  the  integration  finite-difference  and  pre- 
dictor-corrector methods  respectively.  These  charts  show  the  main  subroutines 
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whereas  some  third  level  subroutines  have  been  left  out  for  clarity.  It 
can  be  seen  that  all  of  these  programs  utilize  an  automatic  mesh  genera- 
tion which  is  done  in  MESH  and  POINTS.  The  differences  are  in  the  second 
part  of  the  program,  which  starts  with  the  subroutine  DIFF.  In  the  finite- 
difference  approach,  this  subroutine  is  in  the  Loop  1,  which  is  on  the 
time  interval.  Inside  of  DIFF  there  is  a second  loop  called  Loop  2.  In 
predictor-corrector  method  Loop  2 checks  for  yield  only.  In  the  predictor- 
corrector  method  the  subroutine  DIFF  is  not  in  a time  loop  since  the  sub- 
routine DHPCG  automatically  cycles  through  all  the  time  intervals. 

In  order  to  illustrate  the  actual  program  arrangement,  a listing 
is  given  in  Appendix  C of  the  program  corresponding  to  the  finite-dif- 
ference method.  The  plate  load  is  specified  by  a Subroutine  DISFOR 
which  has  to  be  supplied  by  the  user  of  the  program.  This  subroutine 
converts  the  distributed  load  to  equivalent  concentrated  nodal  forces. 

This  subroutine  is  called  for  each  quadrilateral  element  in  the  loaded 
plane  of  the  plate  and  it  calculates  the  value  of  the  distributed  load  FF(I), 
1=1,4)  for  each  corner  of  the  quadrilateral  and  for  a given  time  para- 
meter TIME.  The  conversion  to  concentrated  forces  is  then  automatic. 

Solution  of  Laminated  Plate 


In  order  to  check  the  accuracy  of  the  dynamic  plate  analysis  it  was 
decided  to  apply  the  finite-difference  approach  to  the  analysis  of  the 
dynamic  response  of  a laminated  plate  which  was  tested  experimentally  at 
the  Ballistic  Research  Laboratories.  The  results  of  this  test  were  reported 
recently?  where. comparison  was  also  made  with  numerical  calculations  per- 
formed with  a hydrodynamic  code  called  HEMP. 

The  laminated  plate  consists  of  three  layers,  two  of  steel  and.  one 
of  aluminum.  The  plate  has  a rectangular  shape  250  x 500  mm,  and  is 
supported  at  the  four  corners  by  resting  on  four  posts.  The  plate  was 
loaded  by  exploding  a circular  explosive  at  the  center  of  the  plate  and 
the  measurements  consisted  of  optical  measurement  of  the  plate  deflection 
on  the  side  away  from  the  explosive. 

In  order  to  model  the  experimental  configuration  by  the  finite-ele- 
ment plate  program,  it  was  necessary  only  to  consider  one  quarter  of  the 
plate  because',  of  symmetry.  Each  layer  of  the  plate  was  represented  by 
eight  rectangular  elements  as  shown  in  Figure  7.  Although  this  is  a re- 
latively coarse  grid,  it  was  considered  to  be  sufficient  in  order  to  de- 
termine whether  or  not  if  the  present  analysis  would  produce  answers  of  v 

the  correct  order  of  magnitude.  The  pressure  load  was  applied  at  the  four 
nodes  closest  to  the  center  of  the  plate  by  converting  the  pressure  to 
time- dependant  concentrated  forces.  These  pressures  were  supplied  by  the  % 

BRL  and  were  numerically  obtained  from  the  HEMP  hydrodynamic  program?. 

The  four  points  at  which  the  loads  were  applied  are  indicated  in  Figure  7. 


^ajerus,  J.  N.,  and  Karpp,  R.  R.,  "Dynamic  Behavior  of  Multi- 
Layered  Plate  Due  to  an  Intense  Impulsive  Load,"  Proceedings  of  the 
Second  International  Conference  on  Mechanical  Behavior  of  Materials," 
Boston,  Mass.,  August,  1976. 
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Since  a relatively  coarse  finite-element  grid  was  used  this  force  distri- 
bution can  only  be  considered  as  a rough  approximation. 

The  numerical  calculations  were  performed  using  a time  interval  of 
0.25ps,  and  the  calculations  were  performed  up  to  60ys.  At  this  time  the 
load  was  almost  equal  to  zero  since  its  decay  is  quite  rapid.  The  deflec- 
tion of  the  plate  was  calculated  by  averaging  the  deflection  at  the  four 
nodal  points  at  which  the  forces  were  applied  and  which  are  indicated  in 
Figure  7.  It  was  considered  that  this  average  would  be  more  representa- 
tive of  the  actual  displacement,  because  of  the  relatively  coarse  finite- 
element  grid.  The  results  of  these  calculations  are  shown  in  Figure  8, 
where  they  are  compared  with  the  experimental  results  and  the  results  cal- 
culated for  the  bottom  plate  by  the  hydrodynamic  code  HEMP.  The  experimental 
results  are  available  only  for  the  bottom  part  of  the  plate,  which  is  op- 
posite from  the  loaded  surface.  The  results  from  the  present  analysis  are 
shown  for  both  the  top  and  the  bottom  surfaces . Considering  the  coarsness 
of  the  trid,  it  is  amazing  that  good  agreement  has  been  obtained  between 
experimental  and  numerical  results 


IV  CONCLUSIONS 

The  method  of  analysis  developed  in  this  investigation  have  been 
successfully  programmed  for  numerical  calculations.  Two  different 
numerical  integration  of  the  dynamic  equations.  The  relative  merits,  of 
these  programs  can  be  stated  at  this  time  in  terms  of  comparing  applica- 
tion to  relatively  simple  problems.  In  these  applications  it  was  found 
that  the  finite-difference  method  was  the  faster  method  although  the 
disadvantage  of  this  approach  is  that  it  is  much  more  prone  to  numerical 
round  off  errors  depending  on  the  choice  of  the  time  interval  size.  The 
predictor- corrector  method  has  a variable  time  interval  which  is  auto- 
matically reduced  according  to  some  specified  error  limit. 


Figure  1.  Coordinates  used  in  the  Analysis  of  Laminated  Plate. 
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Figure  2.  Nodal  Numbering  System 

for  Quadrilateral  Element. 
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Internal  Nodal  Force 


Time 


Figure  4.  Step-Wise  Representation  of  Internal  Nodal  Forces 
in  the  Finite-Difference  Integration  Method 
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Figure  5.  Computer  Flow  Chart  for  Finite  Difference  Method 
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Figure  6.  Computer  Flow  Chart  for  Predictor-Corrector 
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APPENDIX  A 


DERIVATION  OF  DIFFERENCE  EQUATION  OF  MOTION 


This  Appendix  shows  the  derivation  of  Equation  (34) . Consider 
Equations  (33) 


M {A}n+l  " {Vn+l  * {Vn+l 


(A*l) 


[M]  (A) 


n 


{F  } + {F  } 

In  E n 


(A.  2) 


M U}n.!  = <V„.l  - {Fe}„.1 


CAV« 


Applying  the  kinetic  relations,  Equations  (31)  and  (32)  to  the  displacement 
vector  {A}  gives  for  t = (n+l)h 


<A>n.l  - (i>„  * 7 


(A.  4) 


and 


* Wn  * htt)n  * <7  ‘ » 


8n  {A} 


n+1 


(A.  5) 


Similar  expressions  are  obtained  for  t = nh  by  using  Equations  (29)  and  (30) 


tS)n  ' * 7 «A)n-l  * '*>„> 


(A. 6) 


and 


{A>n  . * h a>n_1  * c|  - 6)h2u}n.1 


+ 6h2{A} 


n 


1 2 

Multiplying  Equation  (A. 2)  by  2(j  - 0)h  gives 

2(y  - 6)h2[M](A)n  = 2(j  - B)h2  {Fj}n  + 2(~  - 


»h  {FE}n 


(A.  7) 


(A.  8) 
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2 

Multiplying  Equations  (A.l)  and  (A. 3)  by  0h  results  m 

6h2[M]{A>ntl  = Bh2{Fj}ntl  * 6h2tFE>n*1 

and 

3h2[M]{A>n_1  = $h  + ^FE^n-l 

Adding  Equations  (A. 8),  (A. 9),  and  (A. 10)  gives 

Bh2[M]{A)n^*2(i.-  6)h2[M]{A)n  + 6h2[M]{A}n_j 

- 6h2{FI)nn  * 2(j  - 6)h2{FI)n  * Bh^F^.j 

* 6h2{FBJn.l  * 2 <7  ' W>2{fe,n  * 6h2{FE>n-l 
or 

[M]{h2[g{A}n+1  + 2(i  - 8)(A)n  V BiA>nl]} 

* 6h2[(FI>n.l  * 'I  - 2)(FI}n  * tFI}n-l) 

* eH2[(FE)nn  * - 2)tFE>„  * {FE,n-l] 

Consider  the  bracketed  term  in  Equation  (A. 12),  and  let  it  be 
by  "Term  1,"  then 

Term  1 = 8h2{A}n+1  + 2(1  - 8)h2{A)n+  6h2{A>n ^ 

From  Equation  (A. 5) 

■ <*>„.!  - <4>„  - hti)n  - $ - B)h2{A)n 

Substituting  Equation  (A. 13)  in  Term  1 results  in 

Term  1 = {A>n+J  - {A>n  - h{A>n  - (y  - 8)h2{A>n 

+ 2(1-  8){A}n+  BtA}^ 


(A.  9) 


(A. 10) 


(A. 11) 


(A.  12) 
denoted 


(A. 13) 
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■ <4>n+l  - <4>n  - h(4>„  * <7  - »h2“„ 

+ $h2{A}  . 

n- 1 

.2  .. 

= {A}  . - {A}  - h{A}  + {A} 

n+ 1 n n 2 n 

. yd),  * bk2^  ♦ to),.,  £ 

■ - {4,n  - h(i,n  * “n-l> 

* t-  6h2U)n  - h2(i  - 6)U)n.j) 


From  Equation  (A. 6) 


t-((A>n  * U}n.p  ° h({i)ji  - 


From  Equation  (A. 7) 


- - h2<7  - ««>„-!  - 


{A} 


n 


Substituting  Equation  (A. 15)  and  (A. 16)  in  Equation  (A, 14)  gives 


Term  1 = {A}  . - {A}  - h{A)  + h{A}  - h{A}  . + {A}  n + h{A) 

n+l  n n n n-1  n-1  n 


and  after  cancelling  terms, 

Term  1 = {A},  - 2{A)  + {A}  . 

n+l  n n-l 

Therefore,  upon  substituting  Equation  (A. 17),  Equation  (A. 12) 
be  written  as 

[M]  { U>ntl  - 2U)n  * 

■ 6h2l{Vn,l*  <F-2>  (FE}„  * {FE?n-l 
* lFI>„.l  * <k  ' 2>‘FI>*  * {FI>n-ll 


(A. 14) 

(A. 15) 
(A. 16) 

- {A>n 
(A. 17) 

may 

(A.  18) 
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Rearranging  Equation  (A. 18)  gives 


[M]U>nn  - 6h2{FI}n,1  = 2[M){4}n 

+ (1  - 20)  h2{Fr}n  - [M]{A)n 

* ®h2tPi>Il.1  - 

* h2t6{FE>n.l  * (1  - 3«{FE}n 


which  is  the  same  as  Equation  (34). 
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APPENDIX  B 

COMPUTER  PROGRAM  INPUT 
CARD  DESCRIPTION 
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TITLE  CARD 


Format  (20A4)  Title  (Title  for  particular  case) 


CONTROL  CARD 

Format  (415) 


Columns  1-5 
6-10 
11-15 

16-20 


PRINT  CARD  1-5 


NUMMAT  (Number  of  different  materials ; 6 maximum) 
NUMLA  (Number  of  layers ; 12  maximum) 

NLINC  (Number  of  load  increments  with  time; 

NLINC  >1) 

IPLOT  XPlQt  parameter,  1 if  plot  parameter,  1 
if  plot  required) 

NPRINT  (Number  of  intervals  between  printing) 


INTEGRATION  PARAMETER  CARD 

Format  (F10.5,  E15.7,  I 10) 

Columns  1-10  BET  (0,  acceleration  parameter  or  Newmark's  parameter, 
0 < 0 <_  .25) 

11-25  H (Time-step  size) 

26-35  LINC  (Magnitude  of  load  increment) 


MESH  GENERATION  CONTROL  CARD 


Format  (515) 


Columns  1-5 
6-10 
11-15 
16-20 
21-25 


MAXI  (Maximum  value  of  I in  mesh;  25  maximum) 
MAXJ  (Maximum  value  of  J in  mesh;  100  maximum) 
NSEG  (Number  of  line  segment  cards) 

NBC  (Number  of  boundary  condition  cards) 

NMTL  (Number  of  material  block  cards) 
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LINE  SEGMENT  CARDS. 


The  order  of  line  segment  cards  is  immaterial,  except  when  plots  are 
requested;  in  this  case,  the  line  segment  cards  must  define  the  perimeter 
of  solid  continuously.  The  order  of  line  segment  cards  defining  internal 
straight  lines  is  always  irrelevant. 

i. 

Format  (3(213,  2F8.3),  15) 

I  coordinate  of  1st  point 
J coordinate  of  1st  point 
R coordinate  of  1st  point 
Z coordinate  of  1st  point 
I coordinate  of  2nd  point 
J coordinate  of  2nd  point 
R coordinate  of  2nd  point 
Z coordinate  of  2nd  point 
I coordinate  of  3rd  point 
J coordinate  of  3rd  point 
R coordinate  of  3rd  point 
Z coordinate  of  3rd  point 
Line  segment  type  parameter 

If  the  number  in  column  71  is 

0 Point  (input  only  1st  point) . 

1 straight  line  (input  only  1st  and  2nd  points).  ' 

2 straight  line  as  an  internal  diagonal  (input 

only  1st  and  2nd  points) . 

3 circular  arc  specified  by  1st  and  3rd  points 
at  the  ends  of  the  arc  and  2nd  points  at  the 
mid-point  of  the  arc. 

4 circular  arc  specified  by  1st  and  2nd  points  at 
the  ends  of  the  arc  with  the  coordinates  of  the 
center  of  the  arc  given  as  the  3rd  point  (delete 
I and  J for  3rd  point) . 

5 straight  line  as  a boundary  diagonal  for  which 

I of  1st  point  is  minimum  for  its  row  and/or 

I or  2nd  point  is  minimum  for  its  row  (input 

only  1st  and  2nd  points) . 

6 straight  line  as  a boundary  diagonal  for  which 
I of  1st  point  and/or  2nd  point  is  maximum  for 
its  row  (input  only  1st  and  2nd  points) . 


Columns  1-3 
4-6 
7-14 
15-22 
23-25 
26-28 
29-36 
37-44 
45-47 
48-50 
51-58 
59-66 
67-71 


Note:  In  specifying  a circular  arc,  the  points  are  ordered  such  that  a 

counter-clockwise  direction  about  the  center  is  obtained  upon 
moving  along  the  boundary. 
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BOUNDARY  CONDITION  CARDS 

Each  card  assigns  a boundary  condition  code  to  a block  of  successive 
nodal  points  starting  with  N1  and  ending  with  N2,  inclusive. 

Format  (215,  1 10) 

Columns  1-5  Starting  node  number  N1 
6-10  Ending  node  number  N2 

11-20  Boundary  condition  code 

If  the  number  in  columns  11-20  is; 

0 node  is  not  restrained  (program  assigns 
automatically) 

1 node  is  restrained  in  x direction 

2 node  is  restrained  in  y direction 

3 node  is  restrained  in  z direction 

4 node  is  restrained  in  x and  y directions 

5 node  is  restrained  in  y and  z directions 

6 node  is  restrained  in  z and  x directions 

7 node  is  restrained  in  x,  y,  and  z directions 
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MATERIAL  BLOCK  ASSIGNMENT  CARD 


Each  card  assigns  a material  definition  number  to  a block  of  elements 
defined  by  the  I,  J- coordinates.  One  card  for  each  layer. 

Format  (115,  3F10.0) 


Columns 


1-5  Material  definition  number  (1  through  6) 

6-15  Material  principal  property  inclination  angle 

BETA  in  X-Y  plane 

16-25  Material  principal -property  inclination  angle 

ALPHA  in  N-T  plane 

26-35  Yield  stress-  in  this  material  layer 


PLOT  TITLE  CARD* 

Format  (20A4) 

Columns  1-80  Title  (Title  printed  under  each  plot) 


PLOT  GENERATION  INFORMATION  CARD* 

Format  (2F10.0) 

Columns  1-10  RMAX  (Maximum  x coordinate  of  mesh) 

11-20  ZMAX  (Maximum  y coordinate  of  mesh) 


Note:  Use  only  if  IPLOT  = 1 (plot  required) 
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MATERIAL  PROPERTY  INFORMATION  CARDS 

The  following  group  of  cards  must  be  specified  for  each  material 
(maximum  of  6) . 

a.  MATERIAL  IDENTIFICATION  CARD 

Format  (115,  F10.0) 

Columns  1-5  Material  identification  number 

6-15  Mass  density  of  material  (if  required) 

b.  MATERIAL  PROPERTY  CARDS 

First  Card 
Format  (6F10.0) 

Columns  1-10  Modulus  of  elasticity,  E^ 

11-20  Modulus  of  elasticity,  Eg 

21-30  Modulus  of  elasticity,  ET 

31-40  Poisson's  ratio,  v^g  ' 

41-50  Poisson's  ratio,  v^. 

51-60  Poisson's  ratio,  vgT 

Second  Card 
Format  (3F10.0 

Columns  1-10  Shear  Modulus,  GNg 

11-20  Shear  Modulus,  GgT 

21-30  Shear  Modulus,  G^ 
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LAYER  THICKNESS  CARD 


Format  (12F5.3) 


Columns  1-5  TH(1)  (Thickness. of 'layer  1) 

6-10  ’ TH(2)  (Thickness "of ^ layer  2) 
11-15..  TH(3)  (Thickness  of  layer  3) 

etc>  up  to  TH  (NUMLA).  V 
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APPENDIX  C 


COMPUTER  PROGRAM 
LISTING  FOR  THE 
FINITE  DIFFERENCE  PROGRAM 
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nnnon 


74  / 74  OPT-1 


FTN  4.6  + 42B 


OGRAM  ARMY  JOB 


PROGRAM  ARM Y J OB < INPUT , OUTPUT* TAPE  5- INPUT#  TAPE6-0UTP LIT#  TAP El, 
1TAPE2»TAPE3) 

4—4  -4  -'4.-4 — 4-4-4— 4 — 4 — 4 — 4 — 4 — 4 — 4 — 4—*— 4 — 4 — 4 — * — 4 — 4-  4 — 4 — 4-4-  4-4 — 4 — 4 — 4 — 4—4—4 — 4- 

aaiESC  FINITE  ELEMENT  STRESS  ANALYSIS  OF  AXISYMMETR IC, 

PLANE  STRAIN*  AND  PLANE  STRESS  SOLIOS  WITH  ORTHOTROPIC* 
TEMPERATURE-DEPENDENT  MATERIAL  PROPERTIES 

4—4 — 4 — 4— 4— 4— 4-~4 — 4—  4 — 4 — 4—4 — 4 — 4 — 4 — 4 — 4 — 4—4—4 — 4—4 — 4 — 4-4-4 — 4-  4 — 4 — 4-  4-  4 4 4 
INTEGER  CODE 

COMMQN/BASIC/VOL*NUMNP»NUMEL*NUMLA 
C0MM3N/NPR/NPRINT  . 

CO  MM 9MV-MAT-P^G  ( 1-2 )-.  E-(-9-*-V2  ) * EE-W 

COMMON/ARG/XXX(lO)*YYY(10)*S(24)*XX(3)*YY(3)* 

1CRZ ( 6*  6 ) * XI  (10)*SIG(12)*N*M 
CQMMQN/NPDATA/X(200)*  Y(  200 ) ,C  ODE  ( 200  ) , NPNUMdO,  20) 

GOMMON^E-L-0  A^VA6^4^(l^)^AL^A4-l^>-*-TH-4-},2-)*-I-M-2-DO-*A-T*44^Fft-I-L-(-l-2-) 

C0MM0N/RESULT/D(6*6>*C(6*6)*CNS(6*6 > 

COMMON/TD/IMIN(100)*IMAX(100I,  JMIN(25>*  JM AX ( 25 ) , MAX  I, MAX J* NMTL, NBC 
COMMON / B ASIC  2 /BET»H* TIME *NLAY, IFLAG, IT* NTIME*LINC 

-- — CO  M MO  N/-S-I&MrA-B-SUM^T£UW46-)-»-S-t&IW-<-l-2-  *-50-,-6-),-D-ST-G-(6-UFT6&0) 

COMMON/OISP/DELN1(600)*DELN(600)*DEL(600)*GNM1C600) *GNM2(600) 
COMMON/MASS /A ( 600), B(6OO),CM<0) 

CQMM0N/FIRST/F0(600)»DER(600 1. 

C-OMMQ  N-/-ML-PAX 

COHMON/DELTRI/DELTAI 4*18) 

COMMON/PLYLD/SIGY112 J*DEPS(6) 

4 4 IMENSI0NfcTITLE£20>A  ^ ^ * + ^ + *44444-4  4^44444 

"C  FeaD  AND—WRITE  CONTROL  INFORMATION 

0444444444444444  4-4444444444-4444444444 
50  READ! 5»1000)TITLE»NUMMAT*NUMLA,NLINC* IPLOT 

* &A04-S.*4-OOA-)-N  PRINT 

IF(E0F(5>)920*99 

99  WRITE(6*2000»TITLE*NUMLA*NUMMAT,NLINC 
RE  AD ( 5 * 1002 ) BET*  H*L INC 

WR-I-TE4-6>2-0aa4-B.E-Tl*JH-*X-tNC : 

TIME-0.00 

NT1ME-0 

C4  *4  + 44444444444444444444444444444444 

C &ENER.A-T-E— F-I-NLT-E— EL&MCN-T— HCSM i - 

C4  44444444444444444444444444444444444 

100  CALL  MESH 

IF  (IPLOT. EQ.l)  CALL  MPLOT  * 

M.P.R..I  N.T.-O : 

DO  230  N«1*NUMNP 
IF (MPRINT.NE  .0)  GO  TO  220 
WRITE (6*2003) 

mpju  N-t-3-9 : : — . 

220  MPRIMT-MPRINT-l 

230  WRITE ( 6*  2004 ) N*X(N)*Y(N) 

IT-NUMNP4(NUMLA+1)4B 

-4.4  0 -MP.  R lN.T>  a 

DO  460  M-l.NUMEL 
IF(MPRINT.NE.O)  GO  TO  450 
WRITE (6,2003 > 

MR-R  IN  T - 5-9 : ^ 

450  MPR IMT -MPR1NT-1  ' 

II-IX(N*U 
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non 


ROGRAM  ARMYJOB 


74/74  OPT  -1 


FTN  A .6+428 


7< 


JJ«IX(N/2) 

KK  = I X (N>  3 1 

- L l » I X t N»  4) ; , 

A 60  WR I TE ( 6 * 2009  ) Nz  ( IX  ( N/ I ) z I •!  / 4 ) < 

£-*  + + + <«****  + + + -*  + + ***  + * + ********  + 

C R E AO  AND  WRITE  MATERIAL  PROPERTIES 

...C#-— * — + + — *._+  _* — *....*_*._* — * — * — * — * *_ 

500  CONTINUE 

00  510  M-l/NUMMAT 

READ (5/1004 ) MTYP Ez C RO (MTYPE J ) 

WRI  TE-(-6->-2-0-t0-)— HTY-P-6-* RO  (MTYPE  ) 

READ(  5/1005  ) (E (J/MTYPE) / J-l/9) 

WRITE (6/2011)(E(J/MTYPE)/J" 1/9) 

510  CONTINUE 

ft E A04-5 / -100 6->-(-TM-(  I ) > 1-1/  NUMUU — — 

WRITE (6/ 1007 ) 

1007  FORMAT ("  THICKNESSES") 

WRITE(6/1006  MTH(  I)/ I-l/NUMLA) 

WR  ITE  (4>/-1008-) — - — ' : 

1008  FORMAT ( " YIELD  STRESSES") 

WRI T6( 6/ 1009) (SIGYI I)/ I-l/NUMLA) 

1009  FORMAT ( " ”/ 5 ( 2Xz E 15 . 7 ) > 

00-8 00  - I«T/-NUttLA : — — 

DO  800  J-l/NUMEL 
00  800  K-l/6 
SIGMAU/ J/KJ-0.00 

800  CONTINUE — 

CALL  INIT 

DO  900  NL-l/NLINC 

IF(NL.GT.l)  GO  TO  721 

—00—720 — AL-l/JiUM  LA i 

ALP  HA (N) -ALPHA (N)/ 57.  29 5700 

720  BETA(N)»BETA(N)/57. 295780 

721  CONTINUE 

. roFM_sWFNESTlUTRlx  : 7 ~ 

DO  850  I - 1/ A 

00-8  50  . J -1/1.8 

D EL  TA( 1/ J ) -0. 00 

050  CONTINUE 
CALL  DIFF 

900..  .CONTINUE 

910  GO  TO  50 

1000  FORMAT(20A4/6I5/ F5. 0/515) 

1001  FORMAT(3F10.0) 

10.02.  .FORMAT.!  E10..5./.E13..7-/-11QJ ; 

1 OOA  FORMAT  (I5/F10.0) 

1005  FORMAT(tFlO.O) 

1006  FORMAT  (12F5.3) 

..  2.0.Q.O.—EORM  AT...  (2Hl_z.2.0AAV ,. - 

1 3 3 HO  NUMBER  OF.  LAYERS — -I A/ 

2 33H0  NUM9ER  OF  MATERIALS IA/ 

3 33H0  NUMBER  OF  LOAD  INCREMENTS -IA/) 

2.Q.0.1..F.DRMA-T!  4.1HQ AC  C ELE.RAXIO  N_R.A  R A MEJLE R/-8£T A- FJ.Q.^57 

1A8H0  TIME-STEP  SIZE/H £15.7/ 

2 51  HO  LOAD  INCREMENT/  l INC-- 110/) 
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'03RAM  ASM Y JOS 


74/74 


OPT*  1 


2003  FORMAT 

2004  FORMAT 

2008  FORMAT 

2009  FORMAT 

2010  FORMAT 
21H  » "M AS  S 

2011  FORMAT- ( 


(35H1  NX  Y 

(I5,2F10.4> 

(51H1 EL I J K L 

< I5,4I4,2F13.3) 

< 1H1, "MATERIAL  IDENTIFICATION 
‘ ' DENSITY  «"» E15  «7  ) 


FTN  4.6+428 


) 


77i 


ANGLE BETA  ANGLE-  ALPHA1 : 


NUMBER  •"»I2/ 


11H 
2 1 H 
3 1H 
■ 4-1-H 
51H 
6 1H 
71 H 


, "MODULUS  OF  ELASTICITY-EN  -",E15.7/ 
, "MODULUS  OF  ELASTICITY-ES  »",E15.7/ 
, "MODULUS  OF  ELASTICITY-ET  -",615.7/ 

» " P 0 1 S SON-* A-H-9“NU!IS-^»-,  E W-.  7 / 

, "POISSON  RATIO-NUNT  -">615.7/ 

, "POISSON  RATIO-NUST  -",615.7/ 

...  , "SHEAR  MODJLJS-GNS  »",E15.7/ 

-8-1H-- -,  '-'-SH  6-AR-l!0&ULU-S--GS-T--Jt#  E1-5-.-7J 

91H  ."SHEAR  MODULUS-GTN  »",E15.7/) 

2016  FORMAT  (26H  THE  SYSTEM  CONVERGED 

2017  FORMAT  (33H  THE  SYSTEM  DIO  NOT 

920  STOP 

ENO 


_ IN  1 2, 1 1H 
CONVERGE  IN 


ITERATIONS) 

I 2, 11H  ITERATIONS) 


ERITY  . DETAILS 

I 

I 

I • 


DIAGNOSIS  OF  PROBLEM 

C QN-7R0 L-  VA R I A B-L^-1N-^ OM MON-&R-E  OL)I-V*tr€-NG-E», 

CONTROL  VARIABLE  IN  COMMON  OR  EQUI VALENCED,  OPTIMIZATION  MAY  06 

CONTROL  VARIABLE  IN  COMMON  OR  EQU IV  ALE NCE D,  OPTIMIZATION  MAY  BE 

CONTROL  VARIABLE  IN  COMMON  OR  EQUIV ALENC ED,  OPTIMIZATION  MAY  BE; 


O L IC'  REF  E R E N C E MAP  ( R - 1 ) 


SN  TYPE 
REAL 
-REAL 

RELOCAT ION 
ARRAY  MASS 

APBAY  MA«vS 

14 

ALPHA 

BET- 

REAL 

REAL 

ARRAY 

) 

ELD. 

REAL 
REAL 
REAL 
o E AL 

ARRAY 
ARRAY 
ARRAY 
-ARSIAYL 

ELDATA 
RESULT 
RESULT 
, A 9 G 

0 

2260 

620 

....  . 0 

BSUM 

CM 

CODE 

D 

REAL 

REAL 

INTEGER 

REAL 

ARRAY 

ARRAY 

ARRAY 

SIG| 

MAS- 

NPD: 

RE4 

REAL 

ARRAY 

DISP 

1130 

OELN 

REAL 

ARRAY 

DIS! 

REAL 

ARRAY 

DISP 

0 

DELTA 

REAL 

ARRAY 

DEL 

REAL 

ARRAY 

PLYLO 

1130 

DER 

REAL 

ARRAY 

FIR 

REAL 

ARRAY 

SI  GM 

REAL — 

—ARRAY- 

—MAT 

REAL 

ARRAY 

M AT  P 



7035 



F 

REAL 

ARRAY 

SK 

REAL  ' 

ARRAY 

FIRST 

3410 

SNM  1 

REAL 

ARRAY 

DIS 

REAL 

ARRAY 

DISP 

1 

H 

REAL 

BAS 

INTEGER 

-IF-L  AG 

IKTE^E-R 

BAS 

INTEGER 



144 

IMAX 

INTEGER 

ARRAY 

TD 

INTEGER 

ARRAY 

TO 

13070 

I PLOT 

INTEGER 

43 


GUTINt  ANGLE 


7A/7A 


opt«i 


FTN  A .6  + A28 


77 


C*- 

C 

c + 


SU3RDUTINE  ANGLE  ( R* Z*RC * ZC* ANG) 

INTEGER  CODE  _ . 

. . . C 0 M M 0 N / 3 A S I CV  V 0 L ♦ N UM  N P > N U M E L > NU  ML  A— - - — 

C0MM3N/ARG/XXX(10)*YYY( 10)  * S(2A )*XX  <3)» YY( 3) * 

._ 1 C OM  MO  NV  N P D A T A / X ( ! 0 0 ) } Y I 2 00  ) * COD  E ( 2 0 0 ) , N PNU  M ( 

COMMON /E L DAT  A/ BET A( 12)* ALPHA (12)* TH( 12 ) » I X ( 200* A ) *M ATR U ( 12 ) 

COMMON/TO/IMIN?10O^»I  MAX  U00^»  JMIM25  )»  JMAX (25) * MAX  I* MAX J* NMTL* NBC 

_ -C  0 M MON  /G  A-S-IC2- AB  E-T  *-H  ,-T-PAE  »-Nl-AY-*-IXUA&,-l-T  rNT-IMfrL-fNG— — 

COMMON/SIGM/BSUM. TSUM< 6 )*  S IGMA 1 12 *50* 6) *D SIG( 6) * Ft  300 ) 

COMMON/O ISP/ DELNK 600 )»DELN< 600)* DEL ( 600) *GNM1( 600 ) *GNM2 (600) 
COMMQN/MASS/A  (600)  * B ( 600 ) * CM  < 8 ) • ^ ^ ^ ... 

FIND  ANGLE  OF  INCLINATION  BETWEEN  0 AND  2*PI  ,*„„„„** 

► + ***************  + + ********  + + **** 
PI-3.1A15927 

Dl»  ( Z»ZC-) 

D2»( R-RC ) 


l 


+ * 


IF(ABS(R-RC).GT.l.E-8) 
ANG«PI/2. 

IF.(  Dl.GX.-U-E^a)— iLE-TUJUJ 

ANG«-ANG 

RETURN 

£***************♦ 

X ALL  GW-CIRCLE—' TO-CROSS-AXIS. 

£*******$******** 
100  ANG°ATAN2 ( 01  * D2 ) 

RETURN 

E NO 


OLIC  REFERENCE  NAP  ( R»l> 


GO  TO  100 


+ + ♦*♦♦♦*♦*♦♦♦♦**  + **♦ 


+ * t t t * * ♦ * + « ♦ M A * ^ * 


SN  TYPE 

p P A 1 

RELOCATION 

A 8 B AV  Hi« 

-ALP-HA 

. R-FAt  . 

ARRAY 

^U)i 

REAL 

...  — 

F • P ♦ 

1130 

B 

REAL 

ARRAY 

mas; 

REAL 

BASIC2 

0 

BETA 

REAL 

ARRAY  , 

ELD* 

REAL 

SIGM 

A A 

C 

REAL 

ARRAY  ' 

RESt 

REAL 

INTEGER 

ARR.AY_ 

__  MASS  - 

,110^ 

,-O.NS-  . 

JLEA1 ... 

_ .ARRAY 

RES.l 

ARRAY • 

NPOATA 

6? 

CRZ 

REAL 

ARRAY 

ARG 

REAL 

ARRAY 

RESULT 

2260 

DEL 

REAL 

- ARRAY 

01  Si 

REAL 

ARRAY 

DISP 

0 

DELN1 

REAL 

ARRAY 

0 IS ) 

REAL _ 

REAL 

_ ARRAY 

33 

01  . 

RE-AI 

• 

- • n 

E 

REAL 

ARRAY" 

MATl 

REAL 

ARP  AY 

MATP 

7035 

F 

REAL 

ARRAY 

S IGi 

REAL 
RF  AL 

ARRAY 

DISP 
B.AS  I C? 

4540 

_„j  4. 

GNM2 
...X  F L AG— 

REAL 

INXEGER- 

ARRAY 

DISI 

.aAS_ 

^ ...  n.  .iN  L H W «*-^**  — ' 

INTEGER 

'"TrraT"" 

TO 

0 

IM  IN 

INTEGER 

ARRAY 

TO  ■ 

INTEGER 

BASIC2 

44 

IX 

INTEGER 

ARRAY 

ELD. 

i 

44 


l 


OJTINc  CIRCLE 


SUBROUTINE  C IRC LE ( ANG1, DELPHI , RSTRT, ZSTRT,RC * ZC , It J » 

INTEGERCQDE 

C OMMON/3AS IC/-VOL#  NUMNP*  NUM£  L>-NUML  A — 

COMHQN/MATP/RO(12 >> E (9» 1 2 ) , EE <9 ) 
COMM3N/ARG/XXX(10),YYY(10)»S(2A)»XX(3),YY  ( 3 ) > 

1CRZ (6>6)»XI <10)*SIG(12),N,M 

— COMMON /N PD ATA  / X f 200)  > Y(  200  ) »€30E  (-200)  >NPNUM( 10>  201 

ALPHA(12),TH<  12),IX(200, A)*MATRIL(  12) 
C0MMQN/RESULT/D(6,6)#C (6,6)>CNS (6*6) 

INI 100)>IM4XC 1001, JMIN(25I» JMAX ( 25 ) >M AX  I, MAX J, NMTL » NBC 


CQMM0N/SIGM/BSUM»TSUM(6)*SIGMA(12»50»6)»DSIG(6)>F(60O) 
£SKK2M^0I§I>/DELN1(600,»  DELN(600)»DEL(600)»GNM1(600)#GNM2(600) 
COMMON/MASS/ A (600 ) » B < 600 ) » CM ( 8 ) 

- - DIME  N S-I  ON-  AR-C-1XL>— A0-)-*Ai(4  0»— A 0.). : : 

„ EQUIVALENCE  ( X ( 1 ) , AR J , ( Y < 1 ) t A Z ) 

C****4**********  + ****t*4*t4t4«4t*** 

C PINO  INTERSECTION  OF  LINE  AND  CIRCLE  • NEW  R AND  Z 

.C-*_*.  .*_»-»  ♦ . ♦ » ♦ + **  + **  + *************** 

ANG1-ANG1+DELPHI 

RR«SQRT< (RSTRT-RC) ++2+( ZSTRT-ZC )**2 > 

AR ( I,  J)»RC+.P.R*COS  (ANG1) 

AZ-(  I * J ) «Z£taR.*S-IN.(  ANG-U 

. RETURN 
END 


OLIC  REFERENCE  MAP  -XR.*1). 


SN  TYPE 
REAL 
REAL 

REAL — 

REAL 

REAL 

REAL 

...  . INTEGER 
REAL 
REAL 
REAL 

_ PEAL  ....... 

real 

REAL 

INTEGER 

...  . INTEGER 
INTEGER 
INTEGER 
INTEGER 

_I.N.T.EGER 

INTEGER 

INTEGER 


ARRAY 

.ARRAY. 


ARRAY 

ARRAY 

ARRAY 

ARRAY 

ARRAY 

ARRAY 

ARRAY 

ARRAY 


ARRAY 


RELOCATION 

MASS  1 A ALPHA  REAL 

F.P.  0 AR  REAL 

- HP-D  ATA U 3 0 • -B -RE  AL 

BASIC2  0 BETA  REAL  ' 

SIGM  A*  C REAL 

MASS  110  CNS  REAL 

,NP.DATA . 62—CRZ REAL 

RESULT  2260  DEL  REAL 

OISP  0 DELN1  REAL 

F.P.  7027  DSIG  REAL 

MAT.P 1.70 — EE  REAl 

SIGM.  3A10  GNM1  REAL 

OISP  1 H REAL 

F.P.  A I FLAG  INTEGER 

- -TO 0 ...  IMIN. INTEGER- 

BASIC2  A A IX  INTEGER 

F.P.  3A1  JMAX  INTEGER 

TO  7 LINC  INTEGER 

— — ARG ......  15  O.A.—M  A T R I L I NT  EGER— 

TD  373  MAX J INTEGER 

ARG  375  NBC  INTEGER 


ARRAY 

ARRAY 

ARRAY 

-ARRAY- 

ARRAY 


..ARRAY. 

ARRAY 

ARRAY 


45 


FTN  4.6+428 


SUBROUTINE  OIFF 
INTEGER  CODE 

COMMON /BA  S-I UM  0 U-NUMNB  ,J4UM  E L ^NU  M LA : 

C0MM0N/NATP/R0< 12)#E<9/12), EET9)  # 

COMMON /ARG/ XXX (10 )> YYY( 10 ) >$ ( 2 A ) , XX ( 3 ) # YY ( 3 ) » 
.lCRZC6#6)#XmOI»SIGU2>*N#M  

COMMON/NRO  A-T  A / XT 200  ) # YT  200 ) , COD 6-1200  ) >NPNUM  (1-0,20-) 

COMMON /ELO AT A/ BETA < 12  )/ ALPHA ( 12 >, TH { 12 ) p I X < 200* A ) > MATR I L ( 12 ) 

COMMON/TD/IMIN? 10CH  ilMAXClOO^fJMINtls^fJMAXtZS ) , MA XI , MAX J , NMTL, NBC 

_c  on  Mim/^AS*c-2-/-&£*  

C0MM0N/SIGM/B5UM, T$UM<6>,SIGMA<12,50,6),OSIG<6),F(600) 

COMMON/D  IS  P/DEINK  600 )*  DELN(600),  DEL (600) ,GNM1(600) $ GNM2 (600) 
COMMON/MA$$/A(600)#B(600),CM(8)  ‘ 

-TI  M E-»  T I M 6-.+H-  — — — 

NTIME-NTIME+1 
CALL  STIFF 
CALL  INT 

- CALL-  STRESS 1 

RETURN 

END 


a L I C-  R E RE  R E NC  E-  M A P— ( ft <U- 


* TYPE 
REAL 
REAL 

--REAL 

REAL 

REAL 

REAL 

REAL 

REAL 

REAL 

REAL 

REAL 

INTEGER 
INTEGER 
INTEGER 
~ INTEGER. 
INTEGER 
INTEGER 
INTEGER 
. INTEGER. 
INTEGER 
INTEGER 
INTEGER 

REAL-™. 

REAL 

REAL 


RELOCATION 
ARRAY  MASS 

ARRAY  MASS 

ARRAY—.* EL  DATA- 

ARRAY  RESULT 

ARRAY  RESULT 

ARRAY  ARG 

Aft  (LAY- -D IS-R 


ARRAY 
ARRAY 
ARRAY 
ARRAY 

ARRAY 
ARRAY 
JURA* 


ARRAY 


ARRAY 


DISP 

MATP 

SIGM 

:OISJ> 

BASIC2 

TD 

ELDATA 

_XD 

ARG 

TD 

ARG 

-B  AS  I C.2. 
NPDATA 
BASIC 
BASIC 

-AR-G 

SIGM 

8ASIC2 


14 

0 

0- 

2260 

620 

0 

-.113-0- 

7027 

170 

3410 

-I- 

144 

5 

3 A 1 

J- 

1504 
373 
375 
3 74- 

6 
3 
0 

1 A 0_ 

30 

1 


ALPHA 

BET 

-asim — 
CM 

CODE 

D 

-4VE-LN - 

DS  13  * 
EE 

GNM1 
,K 

IMAX 

IT 

J M AX 

-L INC 

M ATR  I L 
M AX  J 
NBC 

-NMTL  - 
NT  I ME 
NUMLA 
RO 

W - 

TSUM. 


REAL 

INTEGER 

REAL 

— EA  L 

REAL 

REAL 

REAL 

.__.REA1 

INTEGER 

INTEGER 

INTEGER 

INTEGER.. 

INTEGER 

INTEGER 

INTEGER 

INTEGER. 

INTEGER 

INTEGER 

REAL 


ARRAY 


ARRAY 

ARRAY 

ARRAY 

JLRRAT- 

ARRAY 

ARRAY 

ARRAY 

"ar‘ray~ 

ARRAY 

TrrTy" 


ARRAY 

.ARRAY- 

ARRAY 

ARRAY 
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0UT1NE  DI5F0R 


7^*  / 7<t 


opr-i 


FTN  4.6+A28 


77. 


SUBROUTINE  0 IS F OR 

CQMMDN/3ASIC2/aET*H»TlME*NL4Y*IFlAG»IT*NTIME*  LINC 

COMM3N/FOR/FF<  5 U FC-(  5)  - - - 

COMMON/ ARG /X XX ( 10 », YYY(10)*S < 24 ) » XX ( 3 I * YY ( 3 ) * 
1CRZ(6,6)#XI<10),SIG<12>*N,M 
COMMON/NUMB/NITER*HDAT»TO 

- OIMENSION-DISm) — - 

DIMENSION  FF I ( 19 ) 

OATA  (FF1(I)*I»1* 19  1/1.25 531 35 E-l* 7. 198626E-2 *2. 103 87626-2* 
1 3. 6607670E-3»5.117<t375E-3»  3.5  2087B7E-3»2. A256292E-3* 

- 4-* 


- - 1-  1 . 7650A796-3*  W3A-7-2-2A-2€-3»-l,  05-7-37036-3*8 . A7678956 

1 7,0338441E-4*6.0297041E— 4# 5 . 32 36 83 7E-4* 4 . 8614537E-4* 

1 4.3992239E-4*  3 . 93699 41 E-4 * 3 . 4284 337E-4* 2 . 9198733 Er4/ 

11  NITER-19 

TN  ITER  *»-36-.0£--6 

HDAT-2.0E-6 

TO-O.O 

DO  10  I«1*A 

DI $ T ( I ) « S 0-R  T4 XXX- t-I4-**2-+-Y- Y ¥4 1 4**2-) 

10  CONTINUE 

Y3LA$T-.2362*<TIME*1.0E+06*6.  > 

IFiTIME-TNITER ) 15*16.16 

15  ~.  CALL-  - PO I-M-T44-IJ1 6-,-NUM  *4W — 

DO  12  I- 1 * A 

FF(I)»-(FF1 (NUMI+(FF1(NUM+1)-FF1(NUM))*(T IME-TN )/HDAT)*lA,5E6 

12  CONTINUE 

GO-  TO  -IS 

16  DO  17  I»  1,  A 

FF< I >■-* FFI (NITER )+(FFl (NITER J-FF1 (NITER-1 ) >* (T IME-TNITER ) /HDAT ) 

1 +1A.5E6 

I F.(  FF.U.1h.4-T..-0^0.)-FM-M-«X)  

17  CONTINUE 
IS  CONTINUE 

DO  30  ■ I - 1 » A 

■ I F l D IS  T4I-)-,-G-T ,YBL  A-S.TJ. -F-F-I  I.)  «0.0 


30  CONTINUE 
RETURN 
END 


OL  I C R E FE  RE  N C E~  M A F ( R » l'T 
S 

ojl_ : 


R 


SN  TYPE  RELOCATION 

» n a P r p o 


REAL 

0ASIC2 

62 

CRZ 

REAL 

ARRAY  ARG. 

...  . REAL..  . . 

....  ARRAY 

_ 5._ 

FC - 

REAI ., 

__ARRAY FOR. 

REAL 

ARRAY 

FOR"'"" 

117 

F.F1 

REAL 

ARRAY  ! 

real 

BASICZ 

1 

HDAT 

REAL 

NUMB 

INTEGER 

4 

I FLAG 

INTEGER 

BA  SI 

INTEGER 

._  B A5LLC2 

;,Z  ‘ 

-LINC-  — 

INTFGER  -- 

BASI 

INTEGER 

ARG 

15  <r 

N 

INTEGER 

ARG  : 

INTEGER 

NUMB 

3 

N L AY 

INTEGER 

BASI 

i 


) 

i 
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□ UT IN  E IN  I T 


74/74 


OPT  «1 


FTN  4.6+42B 


77  < 


SUBROUTINE  INIT  . -/ 

INTEGER  CODE 

COM  MON /BAS  IC/VOL#  NUMNP-*-NUM£L  »NUML  A = 

CCMMON/MATP/RO  1 1 2 )>  E (9#  12) » EE<9) 

COMMON /ARG/ XXX  (10  ) » Y YY C 1 0 ) # S ( 2 A ) > XX ( 3 ) # YY ( 3 ) # 
1CRZ(6#6)#XI(1'0)»SIG(12)#N#M 

C0MM0N/NPDAT-A/.X.(2OO)»  Y(  200  ) # CODE  (200  )»  NPNUMI  10#  20-) — 

COMMON/E LDATA/B ETA (12) » ALPHA! 12)»TH(12)»IX(200#4)#MATRIL112 ) . 

COMMON/RE  SULT'/D(6#6)#C(6»b.)«CNS(6#6) 

C OMMON/TO/.IM  INI  100  ),  I MAX  (100),  J MINI  25) # JMAX (25) #MAX  I-#  MAX  J>  NMTL  # NBC 

• C OMMDN/ B6,SIC2/-B6-T  » H »-T  I-M 6 #NLAY-#-I  F L AG#-T-T-»NT-IM E #1 1 NC : 

COMMON/S  fGM/BSUM»  TSUMI6 ) » S I GM A ( 12  #<50# 6 ) #DS I G ( 6)  » F ( 600 ) 

COMMON/DI  SP/OELNK  600 ) * DELN(600) # DE L ( 600) # G NM1 ( 600 ) # GNM2 < 600 ) 
COMMON/MAS $ / A ( 600 ) » B ( 600 ) » CM ( 8 ) 

COMMON-AF-llST  A(L0-<600  )-,OERT600) 

C LET  INITIAL  DEFLECTION#  VELOCITY#  AND  STRESS  BE  ZERO 
C LET  INITIAL  FORCE  3E  L INC  •. 

DO  100  I » 1#  I T 

DELN(  U-O-.QO — 

OE R ( I ) = 0«  00 
GNM 1 ( 11=0.00 
F0(I)«0.00  • V 

100--  CONTINUE ...j;;. — 

RETURN 

END  • 


0 L I C RE F E R E NC  E.  -MA  P—LR- 1 ) — — 

S 


SN  TYPE 

RELOCATION 

• 

REAL 

ARRAY 

MASS 

l<t 

ALPHA 

REAL 

ARRAY 

ELD/ 

REAL 

ARRAY 

MASS 

0 

BET 

REAL 

BAS) 

PPAI 

ARRAY 

ELD ATM 

. 0 

BSUM 

FLEA  L 

■ i 'IvL  M ^ " • ’ • • 

PEAL 

ARRAY 

RESULT 

2260 

CM 

REAL 

ARRAY 

MASS 

REAL 

ARRAY 

RESULT 

620 

CODE 

INTEGER 

ARRAY 

. NPD/ 

REAL 

ARRAY 

ARG 

0 

D 

REAL 

ARRAY 

RESL 

REAL.  .... 

ARRAY — . 

...  01  SR- 

1130- 

—D.ELN 

REAL 

. ARRAY— 

04-SF 

REAL 

ARRAY 

DISP 

1130 

DER 

REAL 

ARRAY 

Fl  RS 

REAL 

ARRAY 

SIGN 

14 

E 

REAL 

ARRAY 

MATF 

REAL 

ARRAY 

MATP 

7035 

F 

REAL 

ARRAY 

SIG!“ 

PEAL 

ARRAY 

FIRST 

GNM1 

JLEA-L-  - 

i.AJtftAJt 

REAL 

ARRAY 

OISP 



1 

H 

REAL 

BASJ 

INTEGER 

A 

I FLAG 

INTEGER 

BASI 

INTEGER 

ARRAY 

TO 

0 

I M IN 

INTEGER 

ARRAY 

TD 

. INTEGER 

B A.SI.C2 

A A 

IX  - 

-INTEGER- 

ARRAY- 

E-LD/ 

INTEGER 

ARRAY 

TO 

310 

JMIN 

INTEGER 

ARRAY 

TD 

INTEGER 

BASIC2 

155 

M 

INTEGER 

ARG 

INTEGER 

ARRAY 

ELOATA 

3 72 

MAXI 

INTEGER 

TD 

INTEGER 

TD 

13  A.. 

N 

..IUIFGFR 

ARG- 

’ INTEGER 

TD 

3 

NLAY 

INTEGER 

BASI 

INTEGER 

TD 

1130 

NPNUM 

INTEGER 

ARRAY 

NPD/ 

48 


OUTINE 


I NT 


74/74  OPT  *1 


FTN  4.6+426 


77/ 


SUBROUTINE  INT 
INTEGER  CODE 

...  COMM0N/3ASIC/V43L>NUMN**NUMEL>NUMIA- — 

COMM3N/MATP/RQ(12 ) /E<9»12> # EE (9) 

C0MMQN/ARG/XXX(10)»YYY ( 1 0 > > S < 2 4 ) • XX  ( 3 1 > YY  ! 3 ) > 

1C0MM0N/NPDATA/.X(20o{if(200)»C0D£T200WNPNUM(10»2O). 
CQMM0N/EL9ATA/BETA(12)»ALPHA(12)jTH112)»IX(200>4)>MATRIL(12) 

COMMON /TO/  I MIN ( 100  > * I MAXI 10Q)#JMIN(25>> JMAXJ2 5 > #MAX I#  MAX J;  NMTL#  NBC 
COMMON/BASIC2/BET,H,TIMEjNLAY,IFUG»  lf#NTIME>LINC 

COMMON  / S i GAABSW-f^-S  UMT-6-I » -SlGMA-t-1-2  >5  Ch»  6 > » 0 S IGTfeTj-F^OOT-—-  - — 

COMMON/DISP/OELNl (600)* DELN! 6001# DEL (600)* GNM1( 600) *GNM2 (600) 
C0MM0N/MASS/A(600)»B(600)>CM(8) 

COMMON /F IRS T/FO! 600 )»  DER< 600) 

COMMGN/MASSl/-X/4J-NAM6aO-) 

0 1 MENS  ION  NC00C31 
8 H 9 BE  T ♦H  + ♦ 2 

I H AND  BETA  ARE  INPUT  VARIABLES  STORED  IN  COMMON 

0IF.-H  + *2*{1-2..0<lBE-I-> 

C OBTAIN  INVERSE  OF  MASS  MATRIX#XMINV 
IF! NTIME.GT. 1)S0  TO  201 
DO  200  I ° 1»  I T 

XMI  NVI  I L».l/  A <14 

200  CONTINUE 

201  CONTINUE 
NPT«INUMLA+1)*NUMNP 


I F C NT  I ME-.  GT  . l-J— G0-T0-22-5 ~ , - 

C STARTING  PROCEDURE 
D I F 1 » 0 « 50-B  ET 

C OELN  ANO  DER  ARE  INITIAL  CONDITIONS  ON  DISPLACEMENT 
C • GNM  l...  I S - I N P-UT-  -V  E CT-O-R- 
WR I TE ! 6>  999 > 


799 

.1000 


..206- 


210 


FORMAT!"  INITIAL  FORCE  VECTOR”) 

WRITE(6, 1000) (FO(I)iI»l#IT) 

FORMAT.!  9E 12.  6)  : — 

00  215  K»  1>NPT 
1 0*  COOE ( K ) 

DO  205  J * 1>  3 

-NCOD!  J)*0 — - 

IF(ID.E0.1.0R.IO.E0.4.0R*ID.EQ.7)  NCOOm-l  , 

IFII0.EO.2.QR.ID.EQ.4.0R. ID. EO. 5. OR .ID . EO. 7)  NC0D(2>-1 
IFUD.EQ.  3. OR. 10.  EQ.5.0R  .ID.E0.6.0R.  I0.E0.7)  NCDD(3  )»1 

NP3-3+K.  • 

DO  215  J a 1 » 3 
I-NP3-3+J 

1 F ( NC3D  < J ). E 0 • 0 ) GO  TO  210 
OEL  (I ) aC. 00. 

GO  TO  215 

DEL  II  ) aDE  LN( I ) + H*0  ER ( I ) 


DEL! I ) aOEL  ! I ) + XMINV(I)*!BH*BII)-DIFl*H  + *2+GNMHI ) ) 

DEL!  X.)  “OEL.CI  ).+  t Dl£X*H*.*£.*XMINV(  I ).*F0  !.l ) ) 

215  CONTINUE 
GO  TO  500 
225  00  250  I»liIT 

DELNl.U-)aDEL.NUi : — j — 

DELN ( I ) » DE  L 1 I ) 

250  CONTINUE 
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OUTINE  INT 


270 


74/74 


OPT  *1 


FTN  4.6*428 


77, 


NCOD ( 2 )*1 
NC  OOT-B -)■'!- 


275 

300 

500 


DO  300  K»1,NPT 
I D*C3D £ ( K ) 

DO  270- J -1,3 - ' 

MC0D(J)«0  ", 

I FIID.E3. 1. OR. ID. E 0.4. OR. ID.  EQ.7)  NC0D(1)»1 
IF  ( ID.E0.2.0R. ID.EQ.4.DR. ID. '6 0,5. OR, ID. EQ.7) 

I F ( 1 0 . E 0.3 . OR  .1 D . E 0 .5-.  Oft  .1 0 . E Q . 6 . OR . 1 0 .6  Q . 7 ) 

NP  3s  3*K 
DO  300  J«l»3  V : 

I-NP3-3*  J •• 

- hf  i N G QD  U )-.  64v9->— £0— m-£7-5 : 

DEL  ( I ) *0. 00 
GO  TO  300 

DEL(I)«DELN( I > +2-DELN1 ( I ) 

(>.EL-(  I-)-»0^  ^H^Hi^-(-I4*4a44*a-(-I-)^t4«-l-{-I  )-*D-TF+&H^GHH2-«-H- 

CONTINUE 

CONTINUE 

RETURN  ' ' 

-END , : 


OLIC  REFERENCE  HAP  (R«l) 


SN  TYPE 

K£.AL- - 

REAL 
REAL 
REAL 

I-NT  EGER 

REAL 

1 REAL 

REAL  . 

RIAL- 

REAL 

REAL 

REAL 

- INT-EGER- 

G INTEGER 

INTEGER 
INTEGER 

INTEGER 

INTEGER 

INTEGER 

INTEGER 

INTEGER - 

INTEGER 
INTEGER 
INTEGER 
INTEGER 

A INTEGER 

REAL 


—A  R-R-A-Y- 


RE  LOC  AT  ION 


ARRAY 

ARRAY 

-ARRAY- 

ARRAY 

ARRAY 

ARRAY- 

ARRAY 

ARRAY 

4RRAY 


-MAS-S- 


MASS 

ELDATA 

SIGM 

-N-P-O-A-T-A- 


OISP 

DISP 

-S-IGM- 

HATP 

FIRST 

DISP 


— -14 AL-PWA- 

0 BET 
244  ..BH 
2260. .CM 
-6-2- — — C-R-Z-- 


-ft-EA-L- 


1130  DELN 


1130 

250 


DER 
DIF1 


-1-4 — E~ 
7035  F 
3410  GNM1 
1 H 


ARRAY 

ARRAY 

.-ARRAY- 


ARRAY 


ARRAY 


BASIC2 

• -2-5-2.:- 

144 

Id - 

IMAX 

TD 

5 

IT 

ELDATA 

253 

J 

TUX- 

310-  - 

-JMIN— 

7 

LINC 

ARG 

150  4 

MATRI 

TD 

373 

M AX  J 

ARG 

375. 

mbc 

3 

NL  AY 

TD 

1130 

NPMUM 

254 

NP3 

BAS-IC2.— 

, , ■ 2 

NUMEL 

BASIC 

i 

NUMNP 

MATP 

24 

S 

REAL 
REAL 
REAL 

-REA4- 
REAL 
REAL 
REAL 

-REAL- 

REAL 

REAL 

REAL 

-I  NT-EGER- 


INTEGER 

INTEGER 

INTEGER 

-IH-T-EGE-JL; 

INTEGER 

INTEGER 

INTEGER 

INTEGE-R- 


ARRAY 


ARRAY 

ARRAY 

ARRAY 

ARRAY 

-Aftft-AY- 


BAS] 

MAS 

-ARG 

DISP 

FIR< 

-WA-T-i 


ARRAY  SIGf* 

ARRAY  DISfj 

• BASj 


ARRAY 


-Aftft-4 


TO 
BASI 


INTEGER 

INTEGER 

INTEGER 

-INTEGER- 

INTEGER 

REAL. 


lay TO— 

BASI 

ARRAY  E^DI 

-to- 

BAS] 

ARRAY  NPD4 

Sill 

ARRAY  ARC 


50- 


1 


3 R □ J TINE  INTER 


7 Am 


OPT  «1 


FTN  A.6+A26 


SUBROUTINE  INTER 
INTEGER  CODE 

COMMON  / 3 ASIC  / VQC, NUMN  PrNUM  6 L,  NUML  A - — 

COMM3N/MATP/RO(12),E(9,12) ,EE(9) 

COMMON/ARG/XXXdO),  YYY( 10) , S ( 29 ) , XX ( 3 ) , YY ( 3) , 

1CRZ (6,6) ,XI< 10),SIG(12) ,N,M 

- C D M MO N / 6 L D A T A / BE-T  A ( 12  1,-ALPHA  (12  I , TH(-1Z) , IX  ( 200,  A ),  MATRIl  ( 1-24- 
C3MMDN/RE  SULT/D(6,6),: ( 6, 6 ) »CNS ( 6, 6 > 

CDMMON/TO/IHINdOO)  ‘ ‘ 

CQMMQN/8ASIC2/B£T,H,  unc 
COMMON  / S-IGM-/-&SU*,-f5UM  (-6) 


Dl,«  I 0,01  (INi  l 0, 0)  ' j 

), I MAX (100), J MIN (25), JM AX ( 2 5 > , MAX  I, M AX J, NMTl,  N8 C 
H,TIME,NLAY,IFLAG,IT,NTIME,LINC  ! 

_ SUM  f6-)-,S  I-GM-Ai  12,5  0,6V,  05 IG  ( 6 ) tH  600 I n-\ 

COMMON/DI SP/DELNK 600), DE LN (600 ) , OE L ( 600) , GNM1 ( 600 ) ,GNM2 (600) 

COMM ON /MASS/A(6 00), 3(600), CM(9) 

01  MEN SION  X(7>, Y(7),00(9) 

-OA-TA— Q07-3^-.1259  3 9 18  0S^W»r4*^44  2^94442  7»flVr-,-^2-S, 

1 .6961A0A7802B, .A10A2619231A/ 

X(7>« (XX(1)+XX(2)+XX(3) J/3.0 
Y 1 7 ) ■ ( Y Y (l)  + YY(2)+YY(3))/3.0 

00  100  I-»-W3 

J-I+3 

X( I)»OO(0)*XX (I )+ (1.00-00(8 ) )*X(7) 

X(J  ) = Q0(9)+XX(I )*( 1 .00-00 (9) )*X (7) 

— ¥■  (-I-)-»-00-+-8  )-*-YY  d-1+4-1-,-0©— 0 W-8-)-)-*+-(-7-) : 

100  Y( J) ■0Q(9)*YY(I )+(l. 00-00(9) ) + Y (7 ) 

DO  300  1-1,10 
300  XI(I)»0.00 

— -A  R E A XX41-)*-{-Y-V-<2-)  -Y-YL3KH-X-X-{  2)  ♦ t Y Y-d  1— YAH  144  +**-(  3 ) ♦ (-YY-W 1 

00  AOO  1-1,7 
XI(l)«XId)  + OQ(I) 

X I-<  2->  « X-I-  (-2-)*  00(1)  » X-  (4-4 


xi(3)«xi  oiiooiimm 
XI(<r)-XKA)+OQ(I)*X(I)+Y(I) 
XI  (5  )«XI  (5  ) + 00(I)  *X(I  > + .+  2 
-AOO-X  I 16)  » XI-(  6 )-+  0 CM-I-4-*  Y (-I4-*  * 2-— 
00  500  I- 1, 10 
500  XI(I)»XI(I)*  ARE  A 
RETURN 

E ND-  - — • 


1B0LIC  REFERENCE  ,MAP  (R-l) 

as : 

TER 


:a 

JM 

)E 


SN  TYPE 

REAL 

REAL 

REAL 

REAL 

REAL-. 

INTEGER 
REAL  ■ 


RELOCATION 
...ARRAY. MASS— 


BAS IC  2 
SIGM 

— MASS 


-..  ARft  AY  . 
♦UNOEF 
ARRAY 


RESULT 


. 14  . 

-ALPHA - 

REAL 

__  ^AJIRAY 

-el! 

1130 

B ~ 

REAL 

ARRAY 

MA 

0 

BETA 

REAL 

ARRAY 

EL' 

44 

C 

REAL 

ARRAY  . 

RE 

4X0-.- 

.CNS. 

— RXAL 

array 

—RE 

62 

CRZ 

REAL 

ARRAY 

Aft 

2260 

DEL 

REAL 

ARRAY 

01 

3R0UTINE 


74/ 74 


4.6+420 


SUBROUTINE  LOAD 
INTEGER  CODE 

C O HMD N /•  BASIC  / V OL-,  NUMNP-*  N U M E L*N  UM  L A— 

COMMON/ A RG/ XXX (10)» YYY( 10)*S(24 ) * XX  < 3 ) > YY  t 3)  * 

lCOMMON/NPD4TA/X( loo l # V ( 200  I*  CODE < 200 ) * NPNUM< 10* 20 ) 


C 0 M MO N / E 1 0 A-T A / 8 E-TAM 2 >, ALPHA  U 21# TH(12)  , IX(200*4> ,M  ATR-IL-(-t2I 

COMMON/3 AS  1C 2/3ET*H»TIMe»NLAY» IfLAGt IT*NTIME#l|NC 
COMMON/SIGM/ BSUM, TSUM < 6 1 * SIGMA (12*50* 6>*DSIG< 61 *F< 600) 
C0MMQN/F3R/FF(5)*FC(5»  . 

0 0-  50  I»-l-*4-T — — ; : 

F ( I ) *0 .00 

50  CONTINUE  ■ 

NUMBER  NODES  OF  THE  RECTANGULAR  ELEMENTS 



I1°IX(N*1 ) 

I2  = IX (N*  2)  • 

I3°IX(N*3) 

14  = I X ( N * 4-) : 

: . DESIGNATE  COORDINATES 

xxxm-xui) 

X X X < 2 ) »X  C 12) 

— X XX(34*-X-(-I3I — — : 

XXX(4>»X(I4> 

XXX(5)«l.0/4.0*lXXXm  + XXX(2)+XXX(3>  + XXX<4n- 

YYY (1 ) »Y( II ) ■ 

YY-YI 2- >.«*-(  1-24 -; 

YYY ( 3 ) » Y( 13 ) 

YYY( 4) «Y ( 14  ) 

YYY15 )»1.0/4.0+(YYY(l)+YYY(2)+YYY(3)+YYY(4) ) 

CALL— D-I SF-Q-R— — - : 

IF<FF(1).E0.0.0. AND . FF ( 2 ).EQ.O.O.AND.FF(3).EO»O.O.AND.FF(4).EQ.O.I 
1)  GO  TO  20  „ „ ' 

) FIND  FORCE  AT  CENTER  OF  ELEMENT 

-£F  4-W-B.l.,-0/-4.^0<M-F^-(.ir.L>-F.F-(-24-*-F-F-{-34-+FF-<  4 L>— 7 

: CALL  LOADS  OF  TRIANGLES 

DO  9 I» 1*  5 
FCID-O.  00 

<)_C0  NT- 1 MUE 

CALL  LOT ( 1*2 ) 

CALL  LOT ( 2*  3 ) 

CALL  L QT ( 3 » 4 ) 

C-AL-L — LOT- (-4*1-)- ; 

DO  21  I • 1*4 


FC  < I >»FC ( I ) +FC ( 5) /4.0 
21  CONTINUE 

C-HANGE— XO— GLOBA-L— F-ORC-ES- 

IIl»NUMLA*NUMNP*3+Il+3 

II2=NUMLA*NUMNP*3+I2+3 

II3*NUMLA*NUMNP*3+I3+3 

II 4 * SLU  M LA  + N UMN-P+3.+I-4.+-3-  ... 

F ( 1 1 1 ) 0 F ( III > + F C (1) 
FUI2>*F(  II2)+FC(2) 

F(  1 1 3 » -F ( 1 1 3 )+FC (3  ) 

FI  I I 4 1 « &.(  1 1-4-)  +-  F-C  ( 41  

20  CONTINUE 
RETURN 


-I-N-W^-D  LR&C  T-  ION 
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i ROUTINE 


7A/7A 


A .6+A28 


SUBROUTINE  LQT(II,JJ) 

INTEGER  CODE 

C OMMON/ B AS IC / VO  L»  N UMNR, NUMt  L » NUML  A - - 

COMHaN/ARG/XXX(10),VYY(10),S(2A>,XX(3>,YY(3), 

1 C RZ ( 6, 6) » XIM0)»SIG(12),N,M 
CQMri0N/N?0ATA/X<  200),  Y(  200  ),  CODE  (200),  NPNUMM  0,20)  • 

. COMMON / ELD AT A/3ETA (12), ALPHA  (12  )>TH(12)-»  IX(200>A-)  #,iATRIL-(12-)- 
C0MM0NY3ASIC2/9ET,H,TIMe,NLAY,IFLAG,IT,NTIME,LINC 
CHMMON/SIGM/BSUM, TSUN(6), SIGMA(12,50,6),DSIG(6) ,F(600) 
COMMON/FOR/FF (5), FC(5> 


0 I M ENS I ON-  001  (-3 ) » 002  ( 3 ) ,-DD-3  ( 3 ) , DD  (-3,  3 > , F 6 ( 3 > , F F F-(-3-> — 
DIMENSION  AA(3),9B(3),CC (3) 

THE  TRIANGULAR  DISTRIBUTED  “ 


DESIGNATE  THE  TRIANGULAR  DISTRIBUTED  FORCES 

FFF(1)«FF(II  ) I 

. FFF  < 2 )»FF-(  J J4 

FFF(3)-FF(5) 

XX(1)-XXX(II)  ; 

XX(2)»XXX( JJ) 

. XJC  (-3  )-■  XXX  (-5) ; 

YY ( 1 ) ® YY Y ( 1 1 ) i 

Y Y(  2 ) ® YY Y ( J J ) ! 

YY( 3) “YYY (5  ) i 

...  A A ( I ) » X X-  (.2  M-Y  Y-t-3  )-»X  X4-34-*  Y Y-(-2-> 

AA(2)«XX(3)*YYll)-XX(l)*YY(3) 

AA(3>»XX(l)*YY(2)-XX(2)+YY(l)  ' ! 

BB ( 1 ) ■ YY  j 2 )-YY ( 3) 

B B ( 2 ) * Y-Y  (-3-)~Y¥-(  1 ) . : 

BB(3)«YY(1  )-YY(2)  1 

CC(3)»XX(2)-XX(1)  , 

CC ( 2 ) " XX ( 1 >-XX ( 3)  ! 

— CC.(  I |.iXt(4.WU2  i : -i 

INTEGRATE  XX  AND  YY  I 

CALL  INTER 
DO  12  I®  1, 3 

DD 1 ( I ) - A A_LI ) ±XI-U-».aa.U-)>  X l (-2-)*  c C-(  I-LDLU34 J 

DD2 (I )>AA(I)*XI (2)+BB(I I *XI (5 )+CC ( I )*X I (A) 

DD3(I)»AA(I)  + XI(3)+8B(I)*XI(A)«-CC(I)*XI(6) 

! CONTINUE 

-..0  0...1B  ..IjlIj.3 ... ! 

DO  13  J- 1,3  i 

0D(  I, J)* AA(I)*DD1 (J ) + BB( I >*DD2( J)+CC(I )*DD3(  J ) 

I CONTINUE  • . I 

CALCULATE.  EQ  UI V ELENX.CONCE  N.T.R.AT  ED  ..FORCES ; !■ 

AREA®. 50* (XX (1  ) + (YY(2)-YY( 3) ) +XX ( 2 ) * ( YY (3 >-YY ( 1) )♦ XXC 3 »* ( YY ( 1 ) 

I “YY  (2  ) ) ) i 

. F.£  ( I ) ®lLoiiA...O*J^REA**.2J_*_(  DQ-tl ,1.) ♦-F.F..EJ JJ.tQD.(.l,.2 JAEE^ .(2Jj*DDi-I# JJ-+.F F E 
1(3))  - 

> CONTINUE  •; 

FC(II)«FC(II)+FE(1) 

FC.(  J J >.-F_C  ( JJ.).+  F_EJ2.J . 

FC ( 5 ) » FC ( 5 ) *FE ( 3) 

RETURN 

END  j 


S3 


ROUTINE  MESH 


FTN  4.6+428 


, SUBROUTINE  MESH 
''INTEGER  CODE 

DIMENSION  AR-H-O, — 40  ),  AZ-(  10r  40  ENCODE  (10,-40) ----- < 

COMMON/BASIC/VOL, NUMNP » NUMEL , NUML A | 

COMMON /MATP/R0(12),E(9,12)»Ee(9)  : 

COMMON/ARG/XXX(10),YYY(10),S(24),XX(3)»YY( 3),  ; 

1-CRZ  (6,61,  XII  ia)>  S I GC-12>  > N»  M ----- -•  -rf-.* 

COMMGN/NPDATA/X(2OO)»Y(2OO)»CODE(2O0),NPNUM(lO, 20)  , , . 

COMMON  / ELD  AT  A/ BETA  (12  ),  ALPHA  (12),TH(  12  ) » I X ( 200#  4 ) , M ATR IL  <1 2 ) 

C0MM0N/RESULT/D(6»6),C(6»6)»CNS(6,6)  J 

- C0MK3N/-T  D/EMIN  (40^4  »-4HAX-(-l  004,4  MI-N(25)>J  MA-X-t-2  5 )-r  M A X-I-r  M-A-)Wr  NM-TL-,  AI8  C 

. CQMMGN/3A$IC2/BET,H, TIME,NLAY»IFLAG»IT,NTIME,IINC  , j 

C0MM0N/SIGM/BSUM»T$UM(6)»SIGMA(12,50,6),DSIG(6  ,F<600)  . | 

COMMON/D  I SP/ DELN1 (600), DELN ( 600 ) # DEL ( 600) , GNM1( 600 >,GNM2( 600)  ! 

~ E 0 U I VA L E N C 1 /rtx f 1)  # ARlt? Y("l  >‘?a1"!  j C I X ’(T7 lT,NCOD E)  • 

rC  MESH  CONTROL  INFORMATION-  ' _ ^ • . I 

READ  (5,1000)  MAXI, MAXJ,N$EG, NBC, NMTL 

WRITE  (6,2000)  MAX  I,  MA  X J , N S EG, NBC , NMT  L . ■ I 

£*$*  + ******#$$  * » * * * * + ♦ + * * * * * * * + * * * *’♦  * 

C IN-I-T-IAL4Z  E , 

C + ***********  * * **  + + + + ****  + + **  + •*  * + + **  + 

ISEG--1 

PI-3.1415927 

00  -110— 1-1  ,-40 : — 

00  100  1-1,10 
NCODE ( I , J ) =0 
A R ( I,  J ) » 0 . 

A Z4-I-,  J )-*0-» — 

JM A X ( I ) - 0 
100  JMIN( I )«MAXI , 

IMIN ( J ) »M AXJ 

140-  I M A X-(  J-)  » 0 — 

g****************************  ♦ ♦.**#**« 

C LINE  SEGMENT  CARDS  ~ A ^ 

£*#**♦**  + *****♦*  + *. *****  + * + ♦*. + +******♦, 

450— I SEG«  I SE  6*  1 ; i 

159  IF(ISEG.EQ.NSEG)  GO  TO  400  - . , 

READ( 5,1001)  I1,J1,R1,Z1,I2,J2,R2,Z2,I1,J3,R3,Z3,IPTIQN 
WRITE  (6#  2001  )Ui  Jl#Rl#  Z 1#  12#  R2> Z2# I 3> J 3> R3# Z3#  IPTIDN 

— : — — — 

ARIIl, J1)»R1 


A2CI1fJ1)bZI 

NCODE ( II# Jl)-1  I 

C41-L— *N.I  MXU-T,- J-14 — i 

GO  TO  (150, 200#200, 300, 300,200,200),  IPTION 

g****************************  * * + + ♦*** 

C GENERATE  STRAIGHT  LINES  ON  BOUNDARY 

200  DI-  ABS (FLOAT (12-11) ) j 

DJ»  A3S(FL0AT(J2“J1) ) 

AP.(  12,  J2)»R  2 

- AZ(  I2-,  J24-Z2 - H 

NCODE (12, J2)-l  . 

CALL  MN  IMX  ( I 2, J 2) 
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non  ooo 


(ROUTINE  *ESH 


7A/7A  OPT  »1 


FTN  A. 6 + <»28 


ISTRT-I1 
ISTP« 12 

JSTRT-Jl  - 

JSTPb J2 

0 1 F F«  MAX  1 ( DI  / DJ  ) 

ITER*DI FF-1 . 

IINC«0 ~ - 

J I N C*  0 

IFCI2.NE.il>  1 1 N C » ( 1 2 — 1 1 > /I4BSII2-I1 ) 

IFIJ2.NE.J1)  JINC»( J2-J1)/IABS ( J2-J1) 

1 F ( 1 2 ! N eTiuInoTj 2."n E 7j  1 .TnoT.Tp TI On’.'nI".' 3 f K AP PA «T 

IF(KAPPA.EQ.2)  DI FF»2 . *0  IFF 
RINC- (R2-R1 ) /DIFF 
Z I N C-»-(  Z Z--1 1-  )-/ -D-I-F-F-  - 


WRITE (6/2002 ) DI/ DJ/O IFF /RINC / ZINC/ ITER/ I INC/ JINC /KAPPA 
CHECK  FOR  INPUT  ERROR 


IF (KAPPA.NE. 2.QR.0I .EO.DJ)  GO  TO, 210 
WRITE (6/2003) 

GO  TO  150 

I N TER P 0 L ATE  “ “ 

210  I«I1 

WRITE  (67 200 A j " 

00  230  M «1/ ITER 

IF (ITER. EQ.O. AND. IPTI0N.EQ.2)  GO  TO  230 

I FT  ITER.  EO-.O-^ANO  .LPT-IUN.  EQ-.6J—J&0— TO— 23-0- 

IF( ITER.EQ.O. AND. IPTION.E0.7)  GO  TO  230 

IF( KAPPA.E0.2 ) GO  TO  220 

10LD-I 

.1-1+ 1 INC : 

JDL  D=  J 
J«  J + J INC 

AR(I/J)»AR(IOLD/JOLD)+RINC 

A Z (.  I /- J ) » . AZ-C-IOL-D-/  J.QL  0-1+  Z.I  NC_u — 

WRITE (6/2005)  I / J/ AR { I, J ) , AZ ( I / J ) 

CALL  MNIMXU/J) 

NCODE ( 1/ J ) » 1 

GQ._TCL.230 — 

220  CONTINUE 

IFCI1.GT. I2.AND.IPTI0N. E0.7)  GO  TO  221 
IF(I1.LT.I2.ANQ.IPTI0N.EQ.6)  GO  TO  221 

I ° I +1  INC  ‘ 

AR(  I, J)»AR(IOLD/J >+RINC 
AZ(I/J)“AZ(IOLD/J)+ZINC 

WRIIE..(.6./  200J.) . L/.J./ ARC  1/  JJ./ AZd/J  ) 

NC  ODE ( I / J ) ■ 1 
CALL  MNIMXd/JI 
JOLD-J 

=..  . J - J ♦.  JLIHC : : — 

AR(I. J)«AR(I/JOLD)+RINC 
AZ  ( 1/ J ) ° A Z ( I/JOLDJ  + ZINC 
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BROUT INE 


MESH 


7A/7A  opt  *>i 


FTN  A.6+A23 


NCQDE ( I » J ) * 1 

WRITE (6, 2005  ) I , J » AR ( I, J ) , AZ ( I , J 

C ALL  ■ MNI  M X-M  , J ) - 

60  TO  230 
221  JOLO-J 
J » J + J 

— AR ( I , J ) ® A R-(  I-r  J OLD  ) ♦R-1NC- -- 

AZU»J)«AZ(I,JOLD)+ZINC 
NCQDE ( I > J ) » 1 

WRITE (6,2005  ) I, J , AR ( I , J ) , A Z ( I , J ) 

C ALL-  -MN  IrtX-(-I-rJ-) 

IOLDbI 
I»I  + I INC 

AR( I, J )«AR(IOLO, J ) + RINC 

A Z-(  I r 

NCODE ( I, J ) H 1 

WRITE (6,2005  ) I,J,AR(I,J),AZ(I,J) 
CALL  MNI M X ( I , J ) 

•2-30  CONTINUE- 


IFtKAPPA.EO. 1)  GO  TO  150 
IF(I1.GT.I2.AND.IPTI0N.EQ.7>  GO  TO  231 
IF( Il.LT. 12. AND.IPTI0N.E0.6)  GO  TO  231 

IOLO»-1- 


I * I + 1 1 NC 

AR ( I , J ) ■ AR  ( I OLD , J ) + R INC 
J )+ZINC 


C + 
C— 

c* 


AZ(I, J)»AZ<IOLD,J)*ZIN 

GO— TO-2-32 

231  CONTINUE 
JOLD»J 
Jb J+JINC 

AR-U  ,-J)  b-AR-(-I^-JGLB-)-+RI  NC- 

AZ( I, J)» A Z ( I, JOLO)+ZINC 

232  CONTINUE 
NCQ0E(I,J)b1 

W R I-TE.(  6 , 20054— l^J-,  AR-(-t,-J  ) 

CALL  MNI MX ( I , J ) 

GO  TO  150 

>*<  « ^*<t**#<t%*  + *****  *..+  + **  + * + *******41 


GENE  RATE -CIRCULAR-  ARCS-ON  -BOUNDARY-  - 

******  * + ******************* 
300  AR(I2,J2)»R2 
AZ  t 12, J2)-Z2 
NCODE.Ui>  J 24 — b — 1— 


♦ ♦ 


******** 


CALL  MN IMX ( I 2 , J 2 > 

I F ( IPTION.EQ, 5)  GO  TO  320 


AR ( 13, J3  ) bR3 
AZ(I3, J3)»Z3 
- NC00E-L13.,J  3 )a.l — 


CALL  MNIMX(I3,J3) 
SLAC»(Z2-Z1)/(R2-R1) 
SLBFb-1./SLAC 
S L C E-®..(  Z 3 Z 2 .)./•  ( R3-  R 24- 
SL0F--1. /SLCE 
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OPT-1 


FTN  A.6+A28 


CHECK  FOR  INPUT  ERROR 

IF  ( 49 S ( SI  AC^SICE)  -G-TU  .001 ) -GO-  7 0—310 

WRITE (6,2006)  R1,Z1,R2,  Z2> R3> Z3, SLAC » SLCE 
GO  TO  16  0 

310  R4»Rl+(R2-Rl>/2. 

....  . Z4*ZL+(Z2~Z1 ) /2.  

R5*R2+ (R  3-R  2 ) / 2 • 

Z5=Z2+(Z3-Z2)/2. 

8BF*ZA-SL  BF*RA 

B OF  *-Z  5*  S LD-FJfcft.5 

RC=( B3F-BDF)/ (SLOF-SLBF ) 

ZC«SLBF+RC+BBF 
WRITE <6, 2007)  RC»ZC 

_K  AP.P  A *-l- 


320 

'330 


GO  TO  330 

KAPPA»2 

RC-R3 

Z.C-Z3 

ISTRT-I1 
I$TP*I2 
JSTRT-J1 
J STP*  J2- 


3 AO 


RSTRT-R1 
RSTP-R2 
ZS  TRT*  Z 1 

-ZS  TP*  Z2 

CALL  ANGL  E(RSTRT,ZSTRT, RC»ZC» ANG1 ) 
CALL  ANGLE (RSTP,ZSTP,RC,ZC,ANG2> 

I F C ANG2.LE.ANG1)  ANG2-2. 0*PI +ANG2 


FIND  ANGULAR  INCREMENT 
01*  ABS(FLOATUSTP-ISTRT) ) 

DJ*  A BS-(-FLQ A_T(  J ST-P-J S T.RX) ) , 

I INC=0 

IF(ISTRT.NE.ISTP)  IINC * ( I STP-I STRT ) / 1 A9S ( I ST P-I STRT ) 
.1 F ( J3  TRT-.  NE..J_S-TP-1  - JINC-«UJ ST 8^ J S TR TJ /-LABS^-SJlP-^-IS-TRJl)- 
L AMDA*  1 

IF { IINC .NE.O.AND, JINC.NE.O)  LAM0A*2 
D I F F * MAX  1 ( D I » D J ) 

IT.ER’DIEJLrl.. _.-_i 

I F ( L AMOA . € Q. 2 ) DIFF«2.*0IFF 

DELPHI* { ANG2-ANG1 ) /DI  FF 

WRITE ( 6>  2006 ) ANG1, ANG2# DI FF, DELPHI 


CHECK  FOR  INPUT  ERROR 

[F (LAM DA .NE.2 .OR.DI.EO.DJ)  GO  TO  350 

..WRITE  ( 6,200.3.) u 

GO  TO  150 
350  IO-ISTRT  • 

JO* J STRT 

WRIT  E..L6j^J30AJ. — 

INTERPOLATE 
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9 ROUT  IN E MESH 


7A/7A 


OPT  “1 


FTN  A .6  + 428 


359 


360 


NPT=IAB$(  12-11 .)  +1 ABS  < J2-J1J-1 

DO  380— Msl»  I TER -■*-» 

IF  (LA-1DA.E0.2)  GO  TO  360 
I-IO+IINC 
J-JO+JINC 

CALL  -MNIMX-I  

NCODEU#  J)».l 
CALL  CIRCLE (ANGl# DELPHI *R$TRT# ZSTRT»RC#ZC# I# J) 
WRITE (6»2005 ) I»J»AR(I#J)#AZ(I»J) 

•GO  TO -370 ■ — : 

Is  1 0+1 INC 
J»JO 

NCODE ( I » J ) *1 
-C-ALLMN44XM-»04 


— 3-7-0 
380 


CALL  CIRCLE! ANG1# DELPHI >RSTRT#ZSTRT>RC #ZC# I# J) 
WRITE (6* 2 005  ) I > J # AR  ( I#  J ) # A Z < I # J ) 

J- JO+JINC 

-NCOOE  (I>0-)-=-l 

CALL  NNm<I»J) 

CALL  CIRCLE ( ANG1# DELPHI, RSTRT# ZSTRT# RC# ZC » I # J > 
WRITE (6# 2005 ) I»J#AR(I»J)#AZ(I»J) 

JO  = J 

IFUAHOA.NE.  2)  GO  TO  390 
I*IO+IINC 

-N  G-  QO-E-T-  ) -“-1- 


CALL  MNIMX ( I # J ) 

CALL  CIRCLE* ANG1# DELPHI, RSTRT, ZSTRT. RC# ZC, I* J) 
WRITE  (6*  2 005.)  I , J , AR  ( I,  J > # A Z ( I , J ) 

. 0 —IF-{-K-A  P-R-A-«-i-Q-,-24 — GO — 70—1-50— 

ISTRT-I2 
ISTP= I 3 
JSTRT-J2 
O S TP  «-J-3 — 


399 
C*  * 
-C 


RSTRT-R2 

RST?«R3 

ZSTRT-Z2 

-ZST-P»Z-3 

KAPPA-2 
GO  TO  3 AO 
****** 


******** 


C*  * 
AOO 


******************* 

- C-/LL-C-UL  A-T-E— C-OOR-D-I-NAT-E-S— O-F—IN-LE-RWR — R-O-I-NT-S ——  

******************************* 

IF(HAXJ.LE.2J  GO  TO  A30 
J2=MAXJ-1 

- -DO  - A 2D-N»1  »-50C 

RESIDnO. 

DC1  A 10  J s 2 , J 2 
U-IHINU  )+l 

- 12  = IMAX.LJ  ) -1- 


* * *i 


DO  A10  1“  II# 12 

IF(NCODE(I,J  ).E0.1)  GO  TO  A10  „ . • . ..  14 

DRMARd  + l,  J)  + AR(I-1,  J)+AR(  I,d  + 1)+AR(I,  J-1I)/A.-AR< I# J) 
_0Z*  { AZ-U-+L#J4-+Ai(-T=4»  J-L*AZ-U#-J-+14-+AZ -U-»4»-U.) t#-J4- 
RESID*RESID+ABS(0R)  + A'3S(0Z) 

AR ( I,J)»AR(I, J 1+1.8+DR 
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.ROUTINE  MESH 


7W74 


OPT  “1 


FTN  <i.6+<,28 


A Z ( I > J)«AZ(I, J )+1.8*DZ 

10  CONTINUE  *' 

. - IFtN.EQ.-l)  SES1»RESI0  - - 

IF(N.EQ.I.ANO.RESIO.EQ.O.)GO  TO  A 30 
IF USSID/*ES1.LT.1.E-5)  60  TO  430  . 

420  CONTINUE 

C*^|°*^*^*^*&*^*°  + ,*N*  * ♦ * <i  * * * + **  + **.***♦'***  + **  ♦ 
C ALL  POINTS 

**«**  + ***  + + + + ♦*.  * + *♦**♦**,♦♦***♦♦ 

1000  F 0 R M A-T  • -(  5 1-5 ) : V — ' 

1001  FORMAT  { 3 ( 21 3>  2F8 . 3 )» 15 ) 

2000  FORMAT'  (30H1  MESH  GENERATION  INFORMATION// 

1 41H0  MAXIMUM  VALUE  OF  I IN  THE  MESH ; — 13/ 

- 2 -A1H0 — -M A &I-MUM— VALUE— OF— J — I N -TH E-ME-S «=•-«"»— 3 4 — 

3 A 1 HO  NUMBER  OF  LINE  SEGMENT  CARDS 13/ 

A A1H0  NUMBER  OF  BOUNDARY  CONDITION  CAROS -13/  • ‘ 

5 A1H0  NUMBER  OF  MATERIAL  BLOCK  CARDS 13///): 

2001  FORM  AT- -(-/-/-B 8 H — -INP-U-T 11-  4-1 Ri-  - ---  Ht~,H *7T^TT R2 

12  13  J 3 R3  Z 3 IPTI0N/8X#3C2I4,2F8.4) >16) 

2002  FORMAT  <5H  DI«FA.0»5H  0J»F4.0*7H  DIFF»F4.0*7H  RINC»F6.3»7H 
1NC«F8.3»7H  ITER  « 13#  7H  IINC»I3>7H  JINC»I3,8H  KAPPA-11) 

2003  F-OfLMA  TUT,  38 H * *BA D-I N PU-T— -THI^-L- WE— I-S--NO-T-D-I AG0NA1-) 

200A  FORMAT  (30H  I J AR  AZ ) 

2006  FORMAT  (llH,++^BAD^ INPUT  - THESE  POINTS  DO  NOT  DEFINE  A CIRCLE*/ 

X3X  ^ p | p ^ ^ q y 2^ E 2 0 ______ 

2667  FQRMAT(19H  cente  r c oord  inaTeTT f i V. 6 * i X , F 1 1 . 6 » 1 X ) ) _ , 

2008  FORMAT  ( 7H  ANG1«F9.6,7H  ANG2»F9.6,7H  DIFF-F3.0*9H  DELPHI-F9. 

2009  FORMAT  I//30H  COORDINATES  CALCULATED  AFTER  I3»11H  ITERATIONS) 

_ - RE-TURN — ,~ 

END 


Z 

ZI 


9 

6) 


c VERITY  DETAILS 

i 

I 


DIAGNOSIS  OF  PROBLEM 

CONTR  GL--VARI A3LE-  IN  -COMMON  OR  E QU I VZL  FNCE^—rBPXIMZZ  AZIT1N~41AY~-B 
CONTROL  VARIABLE  IN  COMMON  OR  EO UIV ALE NCEO 9 OPTIMIZATION  MAY  B 
CONTROL  VARIABLE  IN  COMMON  OR  E QUIV AL E NC ED#  OPTIMIZATION  MAY  B 


MBQL  IC REFERENCE ..  MAR 


NTS 

SH 


G1 

F 


SN 


TYPE 

REAL 

PEAL 

REAL:. 

REAL 

PEAL 


RELOCATION 
MASS 


ARRAY 


. AJiRA.Y NRDAIA 

ARRAY  MASS 


1A  ALPHA 
1625  ANG2 

.:_3 10 -AZ 

1614*  BBF 
0 BET 


REAL 
REAL 
— -REAL- 
REAL 
REAL 


ARRAY  EL 

ARRAY — MP 

BA 


BRDUTInIE  mnimx 


74/74 


4 .6+428 


SUBROUTINE  MNIMX(I,J) 

INTEGER  CODE 

C 0 M M ON / B A-$  I C /-V  QW-NUMN  P» NUME L , NUM L A 

^□.NMON/ARG/X  XXUO  ) 1 H Y < 10 1 # 1 < 24  >#  XX  ( 3 ) , YV  ( 3 > , 

-AcOMNON/NRoUim  joor»COOE1200  NPNUMl-0»-2^--—  — 

COMMON /ELDAT A/BET A < 12 )# ALPH A< 12), TH( 121 >IX (200,  4 ), MATRIL < 12 ) 
COMNON/TO/IMINao6)»IMAx1lOO^»  JMIM25  )»  JMAX{25)#MAXI»MAXJ/NHTL>N0C 

COMMON/OISP/OELNK  600) » DELN(60Q)#DEL(600)  # GNMK600  )>GNM2 (600) 

COMMON/MA$$/A<600 )jB(600)#CM<8) 

— I-M  J • -L4WM4-N-*  M-)— : 

IF  < J * GT • J MAX ( I ) ) JM A X ( l ) ° J 
XF(  I«LT»  IMINCJ  ) ) 

I F ( I « GT » I M AX  ( J ).)  IMAX(J)«I 

RrE-TURN 

END 


MBOLIC  REFERENCE  MAP  ( R»  1 ) 


-S4--  TY-R-E 

REAL 

REAL 

REAL 

-- REAL 

REAL 

REAL 

REAL 

fcEJtl 

REAL 

REAL 

REAL 

INTEGER- 

INTEGER 

INTEGER 

INTEGER 

INTEGER 

INTEGER 

INTEGER 

INTEGER 

.INTEGER 

INTEGER 

INTEGER 

INTEGER 

REAL-- 

REAL 

REAL 


ARRAY 

ARRAY 

ARRAY 

_._.-A&*A-Y- 

ARRAY 

ARRAY 

ARRAY 

ARRAY 

ARRAY 

ARRAY 

“ARRAY 


-ARRAY- 


R t LQ-G-A  T I- ON * 


ARRAY 


A.RR.A.Y- 

ARRAY 


MASS 

MASS 

ELDATA 

— -RE-SUL-T 

RESULT 

ARG 

DISP 

0J-5-P- 

MAT  P 
SIGM 
OISP 

— — ■■■—  -F-.-P-  * — — 
TO 

0AS"IC2 

F.  P. 

ARS_ 

TO 

ARG 

— -BAS-IC2 

NPDATA 
BAS  IC 
BASIC 

• ARG 

SIGM 

BASIC2 


1*  ALPHA 
0 BET 
0 BSUM 

-2-260 — CM 

620  CODE 

0 D 

1130  DELN 

-702-7-  -DS I G 

170  EE 
3410  GNM1 

1 H 

A — IHLAG — 

0 IMIN 
44  IX 
341  JM AX 

— 7--  -L  INC-  — 

1 S04  MATRIL 
373  MAXJ 
375  NBC 

374—  NMTl 

6 NTIME 
3 NUMLA 

0 RO 

140 — 5-1 G- 

30  TH 

1 TSUM 


REAL 

REAL 

REAL 

REAL- 

INTEGER 

REAL 

REAL 


IN-WGtft — 

INTEGER 

INTEGER 

INTEGER 

twiiG-ER — 

INTEGER  . 

INTEGER 

INTEGER 


INTEGER 

INTEGER 

REAL 

R-EA1 

REAL 

REAL 


ARRAY 


ARRAY 

ARRAY 

ARRAY 

ARRAY 

ARRAY 

ARRAY 


ARRAY 

ARRAY 

ARRAY 

'array 


ARRAY 
AJRRAY — 
ARRAY 
ARRAY 
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BR3UT INt  MODIFY 


74/ 7^.  OPT  *1 


FTN  4 «6> 423 


!S!iSS5‘?SD?aDIFy 

?K3K«if«a?te?y??5;SKIi.?r'*'  ~ - - 

JAT4/aif??i ?** ApPM5999f t?®0,>NPNUM < 10* 20,  — 

^22Q^«^ULT/Df6f  6>?CC6,6I*CnI(6I6{12,,IX*200,^>'MATRIL<12> 

COMMON/B  ASIG2/B  H-T^  (+J-T-IME-/  N LAY  ,*I- FLAG  ll  T»  NTIME  2L-I N "*X**H*XJ»NHTL,NBC 

COMMON/oisP/DELNifioo{6nF|S-5?Anrt  i2n50,6l>DSIG*6>:'  F(600> 

. C0MM0N/M4SS/M600  ?B?600),LCMm  ' 0E  L ( 600  } > G^Ml  ( 600 ) ,SNM  2 ( 6°0 » 

- C~  A D D ■■  ^ J ^ $ '"'j  j I-T-R-E5-S 

ioo  continue^^ *SIGMA(M'N». n+osica ) 

- CALL  RES L-T 

RETURN  — 

END  i 


MB  OL  IC  R EF  ER  ENC6-- HAP— 


NTS 

DIFY 


TA 

S 

Z 

L 

LN1 


M2 

AX 

AX- 

NC 
TRIE 
X J 

C ..  ... 

TL 
I ME 
MLA 


SN 


TYPE 
REAL 
REAL 

REAL 

REAL 
REAL 
REAL 

REAL- 

REAL 
REAL 
REAL 

REAL - 

INTEGER 

INTEGER 

INTEGER 

INTEGE-tt 

INTEGER 

INTEGER 

INTEGER 

INTEGER 

INTEGER 

INTEGER 

INTEGER 

- REAL  

REAL 

REAL 


ARRAY 

ARRAY 

-ARRA-Y- 

ARRAY 

ARRAY 

ARRAY 

•ARRA-Y- 

ARRAY 

ARRAY 

ARRAY 

AR-ft-AAC 

ARRAY 

■ARfU-Y. 

ARRAY 


ARRAY- 

ARRAY 

ARRAY 


RELOCATION 

MASS 

MASS 

EL-DA-T-A- 

RESULT 

RESULT 

ARG 

OTSP 

DISP 
M AT  P 
SIGM 

-D-IS-P- 


TO 

8ASIC2 
-T-0- 
BASIC2 
ELDATA 
TO 

■TD 

TO 

0ASIC2 
BAS  IC 

-HA  TP- 

ARG 

ELDATA 


1 A ALPHA 
0 BET 

— 0 B-SUH— 

2260  CM 
620  CODE 
0 D 

-a  130 — D*  in- 
7027  DSIG 
170  EE 
3A10  GNM1 

1 — H 

A IFLAG 
0 I MI  N 
AA  IX 

-JMI-N— 

155  M 
372  MAXI 


1130  NPNUM 


2 

1 


NUMEL 
_ NUMNP 
-2  A — S-- - 

7 SIGMA 
2 TIME 


REAL 

REAL 

-.-REAL- — 
REAL 
INTEGER 
REAL 

-REAL 

REAL 
REAL 
REAL 

-REA-L- 


ARRAY 


EL 

BA 

-S-l 


INTEGER 

INTEGER 

INTEGER 

-IN-TE-G-ER- 

INTEGER 

INTEGER 

INTEGER 

-INT-EGE-R- 

INTEGER 

INTEGER 

INTEGER 

-R*AL 

REAL 

REAL 


ARRAY  MA 

ARRAY  NP 

ARRAY  RE 

-AR-R4-Y DI 

ARRAY  SI 

ARRAY  . MA 
ARRAY  DI 

-BA 


BA 

ARRAY  TD 

ARRAY  EL 

-ARRAJf TO 

AR 

TO 

AR 

8 A 

ARRAY  NP 

BA 
BA 

-ARR4* AR 

ARRAY 


ii 
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} 
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HOUTINE  MPLOT 


74/7 4 OPT*! 


FTN  4.6+428 


SUBROUTINE  MPLOT 
INTEGER  CODE 

COMMON/OAS-IC-/-VOL,  NUMNP,  NUMEL»NUML  A r-  - 

COMM0N/MATP/RO(12)»E(9,12),EE(9)  ■ 

COMMON/ARG/XXX( 10) , YYY( 10 ), S (24 j , XX < 3),YY( 3), 

ICR Z (6* 6 )»  XI ( 10 1 • S IG< 12),N:,M 
COMMOM/NPOATA/R (200  ) ,Z< 200) ,C00t<  200 ) »NPNUM< 10,201 


COMMON/Et DAT A/BETA (12)*  ALPHA (12), TH( 12) #1X1200*4 ),MATRIL< 12) 
COM^ON/TD/IMIN?100?>IMAX(100^»lJN^N?25)»JMAX(25)»MAXI#liAXJ»NMTL»NBC 

5§E!  f 1 f o I 

COMMON /D I $ P/ DELN1 < 600 ), D ELN ( 600 ), OEL< 600 ),GNM1( 600) ,GNM2(600) 
COMMDN/MASS/A(600),B(600),CM (8) 

R € A L - X-(-l-00-)-^-Y-t-KI0-)  ,-T-X-(^4^-T-Y-(-24-^T4-T).-€  {-2.04^-Z-HA*- 

READ  (5,1000)  TITLE, RMAX,ZMAX 
CALL  CCP2SY  ( 0. 7, 0 , 2, 0. 2, T ITL E ,0.0, 80  ) 

CALL  CCP1PL  (0.7, 0.7, -3) 

-TX  ( 14  »0 .00 

T Y ( 1 ) *0 . 00 
TX(2)“RM AX/9.0 
T Y ( 2 ) a RM  AX/ 9 . 0 

Z-M  A-X-a-ZM-A-X  ♦J-Y(2)+2-~0 — 

IF  (ZMAX.LT.17.0)  ZMAX«17.0 
00  100  J“ 1, MAX J 
NSTART»IMIN ( J ) 

NST0P-a-lMAX-(-4-) 

N-0 

DO  101  I«NSTART,NSTOP 
N»N  + 1 

NP^NPNUM-(-I-^J-) 

Y{N)-«R(MP) 

101  X(N>aZ(NP) 

CALL  CCP6LN  ( X, Y,N,1 ,TX ,TY ) 

— 1 O-Q-CON-T I NUfs- 

DO  102  I *1, MAXI 
NSTART<<JMIN<  I) 

NS  TOP  ■ J M AX ( I ) 

oo  ids  7»^TarT,_n rfop 

N-N  + l 

NP  » NP  NUM ( I , J ) 

Y..IN-)  a-R-t-NJU 

103  X(N)aZ(NP)  . . 

CA.LL  CCP6LM  (X,Y,N,  1,  TX,  TY) 

102  CONTINUE 

-CA.LL-CC.e-L2.L- tZJtAX-^O.*,-^)— 


1000 


FORMAT 

RETURN 

END 


( 20A4/2F10.0) 


MB OL  i Cf  R. E F E R E N C E MX?  < R «lY 


J 
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FUNCTION  NODE  7W7A  OPT-1  FTN  A.6+428  7, 

FUNCTION  NODE ( I, J ) * 

INTEGER  CODE  i 

COMMON  / 8 AS  IC  / VQL.-NUMNR-.-NUME-L..  NUML  A- — ! 

COMMQN/MATP/RO( 12). E (9.  12) »EE(9)  . , . , i 

C0MMQN/ARG/XXX(10).YYY(I0)»S(2A).XX{3)»YY(3)»  ■ 

1CRZ (6.6) .XI( 10).$IG(12).N,M 

. . - - COMMON/N£DATA/X-(  20O)  #Y(  200)  .CODE  (200). NPNUM(10»20) 

CQMM0N/ELDATA/6ETA(12).ALPHA(12).TH(12)»IX(200»^)»MATRIL(12) 
COMMON/ RESULT/D(6.6)»C(6»6)»CNS(6»6)  ! 

COMMON/TD/IMIN<100)»  I M AX  ( 1 00  ) . J MI  N(  25  ) . JMA  X ( 2 5 ) . M A X I » M AX  J , NMT  L » NBC: 

- - COM  MON  48  AS-IC2  /-8£*^  ,-T-m6^N4.-A-Y-*-I-F L-AG.-I-T-.N TIME>L-I-NC j 

COMMON/SI  GM/8SUM,  TSUMT6)  > SIGMAU2#  50. 6 ).DSIG(  6)  >F  (600) 

COMMON/D  I SP/DELNK 600). DELN( 600). DEL ( 600) . GNM1( 600). GNM2 (600) 
C0MM3N/MASS/A(600  ). 6 ( 600 ) . CM ( 8 ) 

DO  100  : 

NST  AR  T«  IMIN ( J J ) 

NSTOP° IMAX ( J J ) 

D 0 . 1 0 0 • • 1 1 « N ST-ART+.N  S-T-Q  P_ - 

NODE°NODE  + 1 

IF(JJ.EQ.J.AND.II.EQ.I)  RETURN 
100  CONTINUE 

END  7“r  ' ’ ~ 


M B 0 L I C 

REFERENCE 

MAP  ( R * 1 ) 

1 

NT.S 

DE 

. 

( 

SN 

TYPE- 

RF1  HP ATTHN  ' 

: 

real 

ARRAY 

MASS  , 

1A 

ALPHA 

REAL 

ARRAY 

EL 

REAL 

ARRAY 

MASS 

0 

BET 

REAL 

BA 

T A 

REAL 

ARRAY 

ELOATA 

D 

B SUM 

REAL 

SI 

JtfLA! 

ARRAY  

RESULT.. 

1.^2  26  0 - 

„CM 

..  _..RXA4 

AR RAY Mi 

S 

REAL 

ARRAY 

RESULT 

620 

CODE 

INTEGER 

ARRAY 

NPl 

z 

REAL 

ARRAY 

ARG 

0 

D 

REAL 

ARRAY 

RE 

L 

REAL 

ARRAY 

DISP 

1130 

DELN 

REAL 

ARRAY' 

01 

LNl .... 

REAL  

AR5.A.Y 

.DISP 

- 7 02Z. 

- OSXG — 

REA-L 

_A  RRAY St 

REAL 

ARRAY 

MATP 

170 

EE 

REAL 

ARRAY 

NA 

REAL 

ARRAY 

SIGM 

3 A 1 0 

GNM1. 

REAL 

ARRAY 

DI. 

M2 

REAL 

ARRAY 

orsp 

1 

H 

REAL 

BA1 

I NT-tG  E-R  .. 

..  _F.*  P_*  ....  .... 

...,„A~ 

...IF  LAG 

I.N.TEGEH- 

- &A- 

INTEGER 

IAA 

I MAX 

INTEGER 

ARRAY 

H 

IN 

INTEGER 

ARRAY 

TD 

5 

IT 

INTEGER 

INTEGER 

. ARRAY 

FLOATA 

0. 

J 

INTEGER 

F 

INTEGER. 

3 A 1 . 

J M AX  . 

INTEGER- 

...ARRAY 

—ID 

IN 

INTEGER 

ARRAY 

TD 

7 

LINC 

INTEGER 

6A 

INTEGER 

ARG 

I 50A 

MATRIL 

INTEGER 

ARRAY 

EL 

XI 

INTEGER 

TD 

373 

MAX  J 

INTEGER 

TD 

...INTEGER  . 

..ARG  , 

NR  C. 

TNT PC  Eft  — 

-XD: 

AY 

OF 

INTEGER 

INTEGER 

BASIC2 

37A 

1130 

NMTL 

NPNUM 

INTEGER 

INTEGER 

ARRAY 

TD 

np; 

I 

i 

i 
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'ROUTINE  POINT  74/74  DPT»1  FTN  4.6+428 


SUBROUTINE  POINT(TIME»NUM,TN) 
COMMON/NUM9/NITER»HOAT»TO 

AN-  <TI Medial  /HDAT-+1.0  •— — - 

NUM-AN  !-*.■ 

MUM- ( AN+NUM+1 « 5 ) / 2 
TN»TO+(NUM-1)*HOAT 

RETURN-  - - -----  - - 

END 


•BOLIC  REFERENCE  MAP  (R«l) 


ITS  ■ 
INT 

SN 

iER 

't 

..  .T-Y-P-E-- — - 

REAL 

INTEGER 

REAL 

PPA! 

RELOCATION 

NUMB 
F • P * 

NUMB 

■1 

l 

1 

1 

0 

. 0 

HOAT 

NUM 

TN 

REAL 
INTEGER 
. REAL 

NU‘ 

i 

DC  KS 

f'i  C ft  L—  • ■ — - • I--- — •“  ■“  ™' 

LENGTH 

3 

’length 

25B  21 

•I  D • COMMON  LENGTH 

3B  3 
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BRaUTlNE  POINTS 


7A/74  OPT  “1 


FTN  V.6+A28 


SUBROUTINE  POINTS  ' 

INTEGER  CODE 

C DM  M 0 N / 3 A S I G /-VO  L » N UM  N P * N UM  E L * NU  M L A t - - *v 

1COMMON/NPDaU/X(Ioo{»  Y(  200)# CODE i 200  >*NPNUM<  10^ • 

COMMON/EL  DATA/BET  A (12  )»  ALPHA  (12>*TH(12  )*  I X ( 200*  4)*MATRIL(  12  > 
C0MM0N/SESULT/D(6*6)  * C ( 6* 6 ) > CN S ( 6# 6 > ' £• 

CONMQN/TD/IMIN(100)»IMAX<100)» JMIN(25)* JMAX(25)#HAXI#HAXJ»NMTL#NBC 

COM  MON / 3 A S-IG27  BET-rM*  T I ME*  N L AYr-I  FL  AG  # I T » NT-I  ME-*  L I NG- — 

COMMON/S  I GM/BSUM*  TSUM( 6 >*  SIGMA (12*50*6) *D$IG( 6)* F (600) 

COMMON/DI SP/DELN1 (600 )#OELN (600)* DEL (600)* GNM 1(600) *GNM2( 600) 
CQMMON/MASS/A(600)»B(600)»CM(6)  i 

COMMON  APL-YLD-/S  I-GY-(  12-)-*-D  FP-5-(  6-)- 


C* 

-C  ■ 

C* 


DIMENSION  AX(10» 20 ) * A Y ( 10* 20 ) * BL< ANG( 12)*8LKALF(  12) 

E3UIVALE  NCE  ( X ( 1 ) * AX ) # ( Y ( 1 ) » A Y ) 

* ♦ * ; ♦ ***.**************************; 

E ST  A B L ISH-NOO  AL—  P 01 MT-I-NF0IWA44QN 

*****************************  + ***<; 

NEL*0 
NO  DSUM*  0 

DO-  lOO -U-“1*MAXJ 

NSTART«IMIN( J) 

NSTOP* IMAX ( J ) 

DO  100  I “NST  ART*  NS  TOP 

.i.oo-Noasim^NOo-SAm-'H 

NELSUM“0 
JJMAX-MAXJ-1 

DO  110  J J ° 1*  J J MAX  . . 

MS  -TO  P-  M-lN-044J4AX-W4-)-*-I-MAX(^J^-+-l-)-L-4- 

NSTART»MAXO(  IMIN( JJ)* I MINT JJ  + 1) ) 

DO  110  1 I»N$TART  * NSTOP 
110  NELSUM-NELSUM+1- 
NU  MNP-*  NOD-SUM 


NUMEL“NELSUM 
DO  120  J“1»MAXJ 
NSTART»IMIN(J  ) 

- N STOP  “IMA  X-(U4~ 

DO  120  I“N START »NSTOP 
NPNUM ( I * J ) “NODE ( I*  J ) 

NP “NPNUM ( I » J ) 
.X<MP-MAX4~I-*4-) 


1 20  Y ( NP ) » A Y( I*  J ) I 

* * * * * * * * ft******************  * * * * * * * 

READ  AND  ASSIGN  BOUNDARY  CONDITIONS  _ _ ^ ^ ^ ^ ; 

c *******  ♦ <L-JL_*—  *—  *_*_*_.*_J *-U*UJ * — * ■»  * * A * — *_A| 


C* 

c 


c 

C + 


^ * * * * * * * * * * * ♦ * * * * * * * * * * * * * * * 

NUMN“ ( NUMLA+1 )*NUMNP  'i 

00 — 130— I«-1#MUMN - . i 

CODE ( I ) “0 

130  CONTINUE  1 

IF  (N-3C  .EO.O)  GO  TO  210  ! 

...  DO.  200-1  BC0N-*1*  NBC 

READ ( 5 1 1002  )N1*N2> ICN 
DO  200  I = N 1 * N 2 
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non 


FTN  A.6+A28 


i 


CODE  ( I ) » ICN 
200  CONTINUE 

210-  MPRINT-0 — - — : 

DO  230  J»1,MAXJ 
NSTART«IMIN( J) 

NS  TOP  » IM A X ( J ) 

- DO  230  -1  »NST-ART»NSTOP- ~ - - 

NP«NPNUM ( I » J ) 

I F(  MP  R INT • NE  • 0)  GO  TO  220 
WRITE(6,  20001  ' 

IPS  IN-T-*-5<J 

220  MPRINT-MPRINT-1 

230  WRITE (6,2001)  I » J , NP , COD E ( NP ) , X ( NP ) , Y (NP > 

DO  310  I MTL* 1, NUMLA 

-REA  &-(  5 > 1-0 0-04-H  T-L-^aE-T-A-U-M-T-U  , A UP  H A(-I  MT-t)  >-S-I  G-Y-(  1 W-7U) . 

310  MATRIL (IMTL) "MTL 

WR I TE ( 6, 2003  ) < CODE ( I ) , I • l/NUMN  ) 

2003  FORMAT ( 2513) 

ESTABLISH  ELEMENT  INFORMATION 

************************************ 
J J M A X = M A X J - 1 

MIL -I  7 1 

KTL"  1 

DO  AAO  JJ»1,JJMAX 

...  -NS  TOP  ■ M I N 0 ( J-M  AX-(-  M A-X (.J  J-4-14-1- 1 - 

NSTART»MAXO(IMIN( JJ),IMIN< JJ+1)  ) 

DO  AAO  1 1 "NSTAR  T>  NSTOP 
NEL-NEL+1 

-A  2 0-  -I  =N  P N tlM-U  1,-J-J.-)- 

J"I  + 1 

K»NPNUM(  II  + l, JJ+1  ) 

L-K-l 

I x ( M,i  )»I  " " " ; 

I X ( M»  2 ) » J 
I X ( M>  3 ) “K 

IX.(  M,  A )-»l 

A AO  CONTINUE 

IF ( NUMNP  . 6T ,2000)  WRITE (6, 2002  ) 

1000  FORMAT  ( I5»3F10.0) 

1002.  FORMAT! 2.1  5, 110) : -• 

2000  FORMAT  (59H1  I J NP  TYPE  X-ORDINATE  ' Z-ORDINA 

2001  FORMAT  ( 2 1 5 > 16, 1 1 2, FI 3 . 6 , F1A . 6, E26 . 7, E2 A. 7, EZ A .7 > 

2002.  FORMAT.  -(33  H — BAD-  -INPUT-  TOO -MAN  Y NQO  AL—J10 1 NT-SJ 

RETURN 

END 


MBOLI.C  ..RE.F.e.RENC.E.-MAP — (.R*_L)- 
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BROUTINE  QUAD 


7A/74  OPT  «1 


FTN  A.6+<t28 


SUBROUTINE  QUAD 
INTEGER  CODE 

REAL  NUSN»NUT-N,NUTS,NUNS,NUNT,NU$T  

DIMENSION  DUMMY (6,6) » D U M M Y 1 (6,6 ) 

COMMON/ 3 ASIC/VQl,NUMNP,NUMEL»NUMLA 
COMKON/MATP/RO( 12) , E( 9,  12),  EE (9) 

COMMON / ARG/ XXX-(  10  ),  YY YMO)  > S ( 24  > > X X ( 3 ) > YY ( 3)  > 

1CRZ(6»6)>XI(10),SIG(12I»N,M 
COMMON/NPDATA/X(2O0),  Y( 200), CODE ( 200 ) * NPNUM ( 1 0# 20 > 
COMMON/ELOATA/BETA(12)#ALPHA( 12 ) » TH  ( 12  > » I X ( 200,  4 > , M ATR 1 1 { 12  ) 
COMMON/RE-SUL-4701^  61»-C--(-6»M  rCN-ST6^6->- 


10 


C*  * 

C 

C*  + 

100 


- 1-10 

c*  * 
c 

c + * 


COMMON /TO/  IMIN(100)>IMAX(100).»JMIN(25)»JMAX(25), MAXI, MAXJ,NMTL, NBC 
COMMON/3  A SIC  2/BET* H, TIME# NLA Y, I FLAG » IT »NTIME,L INC 

COMM0N/SIGM/BSUM,TSUM  (b  > , S I GMA  ( 12 , 50,  6 I , D S IG  ( 6 ) > F ( 600 ) j 

-C  OM  MO  NV-D-I-S  P/-D  E-tN-l-(-bQ04-»-OE-ti'K  600-l,-0£-L-(60Q>->-G-NMl-(-600->-,GNM2-<600-) 

COMMON/MASS '/A(600),B(600),CM(8)  1 

IF (NTIME.EQ.1I  GO  TO  10 
REAO(l)  ( ( C R Z C I > J>, J»l,6), 1 = 1,6  ) 

. GO  -TO  -151 
CONTINUE 
I1«IX(N,1 > 

J1*IX(N,2> 

— K-l-a  I-X-(J4>0 1 
L1»IX(N,4) 

MTY?E«MATRIL(NIAY)  1 

**********************  * * ********** 

- INI  E R P OL- A-TE-M  A-1E-R-I Al— P-R-OP  ER-T-l&S— — — ! 

********************************** 
00  100  1-1,9  1 

E £ ( I ) =E  ( I , MTY PE ) 

00— 110-I--1..6 : : 

DO  110  J » 1 > 6 

CNS ( I * J >*0,00  

C(I,J)«0.00 

-D-(-  0 .00 

* * * * * **************  * * ************* 


FORM  STRESS-STRAIN  RELATIONSHIP  IN  N-S-T  SYSTEM 
***********************.*********** 

••  NUNSa-EE-t  44 4 

NUNT-EE ( 5 ) ! 

NUST  = E E ( 6 ) ' ! 

NUSN»(EE (2)+NUNS)/EE(l>  . . i 

... NUTN»4-EE-(-34-*-NUHTU7-EE.m : 1 

NUTS* (EB<3)«N0ST>  i 

0 1 V*l. 00-NUN S*NUSN-NU$T*NUTS-NUNT*NUTN-NUSN*NUNT*NUTS 
i-NUNS+NUTN+NUST  . • 

-C-N-S-l  EiUUAt-L..aO=NUST-+-NUTS.)/4)I V- : J 

CNS(1,2>»EE(2)*(NUNS+NUNT«‘NUT$)/DIV  • 

CNS(1»3)=EE(3)*(NUNT+NUNS*NUST)/DIV  . ‘ 

CNS <2,1 ) «CNS (1, 2) 

- CNS  U,2)  = EE-(-2.»-(  1.00-NUNT*.NUTN  >/  D I V 

CNS  ( 2,  3 ) * EE ( 3 ) + (NUST-*-N'JSN*NUNT)  /OTV 
CNS(3,1)=CNS<1,3) 

CNS(3,2)=CNS(2,3) 

— -C NS  ( 3,  3)  =-££-(  3)-M  1 , 00--r.NU  NS-*NUSFU  /04V : 

CNS(4,4)»EE (7) 

CNS (5, 5>=EE(8) 
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^ROUTINE 


QUAD 


74/74  OPT  *1 


FTN  6+428 


SET  UP6SrRAIN  TRANSFORM  TO  N-S-T  SYSTEM 
. S I N 4«  S I N<  A L-PHAT  M4-I-— 

COSA»COS( ALPHACM) ) 

S2»SINA**2 
C2“C0S A++2 
■ SC“S  INA-+COSA 
D(  If  1 )>C2 
D(  If  3) ”S2 
D(  If 6)»-SC  . 

~D  ( 2-f  14*  S 2 

D ( 2f  3 ) °C  2 
D(2f 6)»SC 
D ( 3 f 2 ) a 1 « 00 
D<  4f-l)  »2-*-0Q*-SC- 


0 ( 4»  3 ) “-2 • OO+SC 
D(4,6)*C2-S2 
D ( 5f  4 ) *S INA 
0(5f5)»CnS-A 


D(6f4>»C0SA 

SEt'u^STRAIN  TRANSFORMATION  TO  R-Z-T  SYSTEM 

S I N3J.  S I iu  B.E.TJU  H44 

COSB-COS (BETA(M) ) 

S2-SINB++2 

C2-C3SB++2 

SC-»Sim*COSB 

ccif 1>«S2  . 

C ( 1»  2)*C  2 
C ( 1.4)«SC 

C (2*.l  )»C-2 

C ( 2f  2 ) = S2 
C ( 2 f 4 ) »- SC 
C ( 3f  3 > = 1 • 00 

— C C.4.»  1 ).«^2V00iLSC 

C (4f2)»2. OO+SC 
C (4f 4J-S2-C2 
C ( 5 f 5 ) “ S INB 

_.C  ( 5 .#  6.)..».=,-C0SR..~  - 

C ( 6/5  ) *C  OSB 
C(6f6)*SINB 

: CALCULATE  CRZ  MATRIX 

..00  12.0_JL“1.  6 

DO  120  J ® 1 f 6 
DUMMY ( If  J )■ 0. 00 
DO  120  K »1 f 6 

_ 1.2  0.. JM1M  M Y.  ( X t-  4.)  -S-DiiMM-Y-C  Lt- J.)  -+  0 •(  If  KJACTKf.JJ- - 

DD  130  1 3 1 f 6 
DO  130  J 3 1 f 6 
DUMMYUIf  J ) “0.00 

130  DUMMY1  ( I * J ) “DUMMY  i C" If  J f+C  N STlVk)  + OUMMY  (Kf  J ) 

DO  140  I “ 1 f 6 
DO  140  J“ If 6 

D.UMM.Y.LIf  JJ_”£L».D.O 

140  DUMMYtlf J)“DUMMY(I» J)+D(Kf I)*DUMMYl(Kf J) 
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BfiUUTINE  QUAD 


74/74 


OP  T*  1 


FTN  4.6+428 


DO  150  1*1,6 
DD  150  J«l,6 

C R Z ( IVJ  ) * &.00- - - - 

00  150  K = 1 *6 

150  CRZ(I»J)=CRZ(I»J)+C(K,I )+ DUMMY (K»  J ) 
WRITE (1)  ( (CRZ< I, J)*  J*1,6),I*1,6> 

151  CONTINUE - ' -••••- 

DO  200  1-1,4 
MM*  I X (N, I ) 

-XXX  I I )-  X-IN-M-) — 

XXXU  + 4)-  X(ilM) 

Y Y Y( I ) ■ Y(MM) 

200  YYY (1+4)*  Y ( MM) 

XXX  ( 9-) «— (X-X-X-(-l-)-+- 

YY Y ( 9 ) * ( YYY  (1 ) + 

XXX(IO)*  XXX ( 9) 

YYY(IO)*  YYY (9  ) 

DO  2 50  -I°l,8 

CM ( I ) -0. 00 
250  CONTINUE 

DO  900  IA*1, 24 

900  S ( I-A-)  *— 0 .00 

CALL  TRISTF ( 1, 2) 

CALL  TRI ST F ( 2, 3 ) 

CALL  TR ISTF ( 3,4 ) 


XXX  (-2-J-+-X  XX-L3-H—  -X  X X-<  4 -B  A4-.  0 Qh 
YY Y ( 2 ) + Y YY ( 3 ) + YYY(4))/4.00 


C A L-L--  TR-I  S-T  F ( 4,1-)- 
RETURN 
END 


MB  OL-I  C-  -R  E F-ER-SNCE  -M  AP-C-R-4-)— 

NTS 

AD 


SN  TYPE 

RELOCATION 

REAL 

ARRAY 

MASS 

REAL 

ARRAY 

MASS 

IA 

fie  a i 

. AftHA-Y- 

-E-LOAXA- 

real 

ARRAY 

RESULT 

S 

REAL 

ARRAY.. 

RESULT 

SA 

REAL 

2 

REAL 

4^^AY - 

ARG 

REAL 

ARRAY 

RESULT 

LN  . 

REAL 

ARRAY 

DISP 

V 

REAL 

MMY — 

REAL  - - 

. _ array 

. 

REAL 

ARRAY 

M AT  P 

REAL 

ARRAY 

SI  GM 

M2 

REAL 

ARRAY 

DISP 

INTEGER- 

. ... 

! 

LAG  ' 

INTEGER 

BAS IC2 

IN 

INTEGER 

ARRAY 

TD 

......  . 

* 



14  ALPHA 


REAL 


ARRAY 


EL 

1A 


0 — 

2260 
620 
515 
■5  1? 

B-SUM — 

CM 

CDOE 
CQSB. 
r ? 

.--R-EAX 

REAL 

INTEGER 

REAL 

PEAL 

ARRAY 

ARRAY 

MA 

NP 

2260 

V U'  — — '■ 

DEL 

REAL 

ARRAY 

D I 

0 

DELN1 

REAL 

ARRAY 

01 

7027 

DSIG 

REAL 

ARRAY 

SI 

565- - 

DUMMY1 — 

....REAL- 

-ARRA.Y — 

..  — ! 

170 

E E 

REAL 

ARRAY 

MA 

3410 

GNM1 

REAL 

ARRAY 

D I 

1 

H 

REAL 

BA 

T A 

.INTEGER 

- • • U — - 

144 

IMAX 

INTEGER 

ARRAY 

TD1- 

5 

IT 

INTEGER 

BA 
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BROUTINE  RESLT 


7A/7* 


OPT  “1 


A .6+A28 


SUBROUTINE  RESLT 

INTEGER  CODE  _ * ■ 

C0MM0N/8ASIC/VQU*NUMNP#  NUMELiNUMLA  - - — 

COMMON/ MATP/RQ  ( 12  ) # E ( 9 # 1 2 ) # E E ( 9 ) ^ k 

COMMON/ARG/XXXUO)#  YYY(10)#S(2A)#  XX(3)#  YY<3># 

XCRZ <6# 6)#XI CIO) #$IG( 12) # N#M 

COMMON/NPDATA/X(  200  ) # Y ( 200  ) t CODE  < 200  ) fNPNUM  ( 

COMMON /EL  DAT  A/ BE TA< 12 )# ALPHA (12) # TH (12 ) # I X < 200# A ) # MATR I L ( 12) 
COMMON/ RESULT/D (6# 6) #C(6# 6 ) #CN$ <6*6 ) 

COMMON/ TO / IM IN (100)# IMA X( 1001# JMIN(25 ># JMAX<25) #MAXI# MA X J # NMTL • NBC 

COMKQ&V£A-SYC4t/4ME-T  *-H#  YIALE-#  NLAY^I  FL  AG  T ^NT-I^ 

COMM0N/SIGM/BSUM#TSUM(6)#SIGMA(12#50#6)#DSXG(6)#F(600)  , . 

COMMON/DISP/DELNK  600 ) # DE L N ( 600 ) # DE L ( 600) # GNM1 { 600 ) # GNM2 ( 600 ) 
COMMON/ MAS$/A(600)fB(600)#CM(8) 

- RE  T URN — H 

END  i 


MBOLIC  REFERENCE  MAP  (R=l) 


SM~  TYP-E 

REAL 

REAL 

REAL 

REAL — t - 

REAL 

REAL 

REAL 

. REAL- 

REAL 

REAL 

REAL 

LMT EGER- 

INTEGER 

INTEGER 

INTEGER 

INTEGER 

I NTEGER 
INTEGER 
INTEGER 

..  INTEGER- 

INTEGER 

INTEGER 

REAL 

REAL  . ... 

REAL 

REAL 

REAL 

REAL  

REAL 


RE-iaCA-T-ION 

ARRAY  MASS 

ARRAY  MASS 

ARRAY  ELDATA 

AR.RA.Y~ JtESULY — 

ARRAY  RESULT 

array  arg 

ARRAY  DISP 

ARRAY-— — CL  ISP : 

ARRAY  M ATP 

ARRAY  SIGM 

ARRAY  DISP 

3 AS-IC-2 

ARRAY  TD 

ARRAY  ELDATA 

ARRAY  TD 

; ARG 

TD 

• ARG 
BASIC  2 

.ARRAY NP.DAT.A 

BASIC 


ARRAY 
, ARRAY. 


ARRAY 

ARRAY., 

ARRAY 


BASIC 
■ ARG 

SIGM 

3ASIC2 

BASIC 

ARG 

ARG_: 

ARG 


1A 

0 

0 

620 

0 

1130 
.7027 — 
170 
3 A10 
1. 

— 1-AA- 
5 

3 A 1 
7 

15  Q A 

373 

375 

37A 

_6.- 

3 

0 

1 AO 

..jo,.: 

l 

o 

5 A 

......310 

12 


ALPHA 

BET 

3SUM 

-C-M 

CODE 

D 

OELN 

— D-SYG r 

EE 

GNM1 

H 

~I  MAX — 

IT 

JM  AX 
LINC 

MA.T.RI.L— 

M AX  J 

NBC 

NMTL 

.N.TJ&E 

NUMLA 
RO 
S 13 

: TH 

TSUM 

X 

XX 

_..Y 1_ 

YYY 


ARRAY 


-R£AL ARRAY — 

INTEGER  ARRAY 

REAL  ARRAY 

REAL  ARRAY 

-REAL ARRAY- 

REAL  ARRAY 

REAL  ARRAY 

REAL 

-INXE  -G-ER ARRAY- 

INTEGER 

INTEGER  ARRAY  . 

INTEGER 

-INTEGER ARRAY — 

INTEGER 

INTEGER 

INTEGER 

, INTEGER..- 

INTEGER 

REAL  ARRAY 

REAL  ARRAY 

. REAL ARRAY- 

REAL  ARRAY 

REAL  ARRAY 

REAL  ARRAY 

-REAL ARRAY — 

REAL  ARRAY 


BROUTINE  STIFF  74/ 7A  OPT«l  FTN  A. 6+428  7 


SUBROUTINE  STIFF 
INTEGER  CODE 

COMMON/BASIC/VOL,-NUMNP,NyME.L,NUMLA  - — 

CQMMaN/MATP/R0(12)»E(9,  12),  EE  (9) 

COMMON/ ARG/ XXX  (10) ,YYY(  10 ) , S ( 2 A ) , XX  ( 3 ) , Y Y ( 3 ) , 
1CRZ(6,6),XI(10),SIG(12),N,M 

COMMON/NP  DAT  A/X(  200)  »Y<  200)  ,COOE ( 200) ,NPNUM  (10,70)  

COMMON/ELD AT A/9  ETA (12 ), ALPH A ( 12 ) , TH ( 12 ) , IX ( 200, A > , M ATR  IL  ( 1 2 ) 
C0MM0N/RE$ULT/D(6»(i),Cl6»6)»CNS(6,6) 

COMMON / TO /IM  IN  (100),IMAX(  100),  JMIN(25),  JM AX ( 25 ) , MAX  I, MAX  J,NMTL,NBC 

• C Oil  MO  N / 8 A S-I  C 2- / B74-,  H »-TIME,NL-A-Y*-I  F L AG, -I  T , NT  I ME  ,1-1 NG 

COMMON/S IGM/BSUM,TSUM(6) ,SIGMA(12,  50,6)>D$IG(6),F(6CO) 

COMMON /D I SP/0ELN1 (600), DELN( 600), DEL (600), GNM 1(600 ),GNM2( 600) 
C0MM0N/MASS/A(600),9(600),CM<6) 

- --C0MM0N./MASSl/-X-M-I-N-V-(-6OO)- — : 

I FL  AG« 1 

IF (NTIME. EO.l ) GO  TO  75 
DO  50  I “1 , IT 

• • - — GNM2  ( I ) ’GNMl-(-l) 

GNM1 ( I ) ®B  ( I ) 

50  CONTINUE 
75  CONTINUE 

- 00  - 100  -Is.1-,4-7 

D(I)«  0.00 

100  CONTINUE 
IF(NTIME.GT,1)G0  TO  102 

DO  101-1-1-,-IT — - — 

A ( I ) - 0.00 

101  CONTINUE 

102  CONTINUE 

REWIND I.—  ..- : 

REWIND  2 
REWIND  3 

DO  3 AO  H» 1, NUML  A 

00  3 AO  N*~ L , NUM E L : * 

CALL  QUAD 
00  3A0  I « 1 , A 

II  . » 3*I.X..(N,  I.)  + 3A(Msl4  «NUMNJ»-  .-: r 

J»  1 1-2 

DO  3 A 0 K = J , II 

J J *K- 1 1 + 3*1 

B (.<)»..  3 ( X > + S (JJJ ' _ 

Ki<«  K + 3 + NUMNP 
B ( KK  ) » 3(KK)+S( J J + 12) 

IF(NTIME.GT.1)GQ  TO  3A0 

. . .AIK  ) ■ A ( ILhC  M.CXJ , 

A(KK)»A(KK)  + CM(H-A) 

3 AO  CONTINUE 

CALL  LOAD  . 

. . DO  A.OO  I=.l  ,.IT _ 

A 00  9< I » * -3 ( I > + F( I) 

RETURN 

END 
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(ROUTINE  STRESS 


7A/7A 


A.6+A20 


SUBROUTINE  STRESS 
INTEGER  CODE 
C OM  MON/  ? A S-I-C-/-VQL-»-NUMNPj 
COMMON/N^R/NPRINT 
COMMON/MATP/ RQ112  ) # E ( 9j 


NUHEL'»NUMIA 


COMMON /M ATP/ RQ( 12 ) # E ( 9# 12 ) # EE ( 9 ) 

COMMQN/ARG/XXX( 10 ) * YYY (10 ) , S ( 2 A ) * XX ( 3 I > YY ( 3 ) * 

1 C R Z < 6 » 6 ) » X M l 0 > , S I G ( 1 2 ) ■*  N , M - - - — — -- 

COMMON/NPDATA/X{20Q)*Y(200)»COD£(200)*NPNyM(10*20)  , 

CONNON/ELDATA/BETAI 12  >»  AL PH A 112 ) » TH { 12  ) » I X ( 200* 4 I # HATR U ( 12) 
C0MM3N/RESUIT/D(6*6>»C(6*6)*CNS<6*6).  . 

COM  MON  / T 0 / IN  I N ( 1 0G0-#I  N AX-(i 0 0>-»  J MIN- (-2-5  )»  J M AX-TS^DM-A-X-I-r  N A'XtL*-NHT'L-/N&C' 

COMMDN/9ASIC2/BET*H*TIME»NLAY#IFLAG#IT#  NTIME*LINC 
C0MMDN/SIGM/9SUM. TSUM  (6  ) * S ISMA  ( 12*  50*6  ) *D$I  G(  6)  # F (-600)  . 

COMMON/O ISP /DEINK  600)*  DELNI 600  I » DE l t 600) *GNM1( 600 ) #GNM2 (600 ) 

COM  M10N-/M-A-SS-/-A  ( ^<W-*-B-<-600  >-*-GM  (-B-) 

COMMON/0cLT/XOEL(A#18) 

COMMON /D  ELTR I /DELTA ( A*  18  ) 

COMMON/PL YLD/ SI  GY ( 12 )*DEPS (6) 

DIMENSIONED!)  E-L  (-6  0 0 ) r — ; : 

IF(NTIME.EQ.l)  GO  TO  96  _ 

I F ( (NTIME/NPRINT)*NPRINT.NE.NTIME ) GO  TO  99 
93  CONTINUE 

WRI-TE-(6>L00-J > N T-I  ME 

1003  FORMATC  DEL  FROM  SUBROUTINE  S TRE  SS*  C YC  LE  ■"  * 15  ) 
WRITE(6,1001)(DEL(I)*I»1*IT) 

1001  FORMAT ( 9E 1 2 • 6 ) . . 

9 9 -C  ONT I NUE : — 

C CALCULATE  CHANGES  IN  DISPLACEMENTS 
DO  100  1*1# IT 
DDEL(I)*OEL(I)-DELN(I) 

. ._4  00-  XON-T-IN  U6 -• 

DO  101  I » 1 » A 
00  101  J * 1 * 18 
X DE  L ( I * J ) * 0. 00 

D£.L-TA-(4->J.)  *0-»00 — ^ 

101  CONTINUE 
IFLAG»2 
REWIND  1 

-R  E W-INO-2 

REWIND  3 

DO  200  M*l#  NUML A 
NLAY-M 

DO  _2QjO.JU.UNU  M E-L : 

DO  150  I • 1 > A 
11=1+1 

I F ( I. ED. A)  11*1 

MM*-IX(-N»4-) ' 

MM1*IX(N>  ID 

11*3+ MM+3* (M-l ) *NUMNP 

II1-3+MM1+3* (M-1)*NUMNP 

J 1*1.1 1-2 ; 7 

DO  1A0  K« J> 1 1 
IK=K- J+l 

X QE  L ( L*  IK-)  aO-DE-L  IK ) : 

OcLTA(IjIK)*DELN(K) 

KK=K+3*NUMNP  ' 
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ROUT  INE  STRESS  74/74  OPT-l  F TN  4.6+420  7 


XDSL!IjIK+6)-D0£L!KK) 

DELTA!  I,IK+6)»0ELN<K!<> 

140  CONTINUE  - ' - - 

DO  145  <»J1,II1 
IK-K-J 1+1 

XDEL!I»IK+3) -DDEL IK) 

-0ELTA(I,IK+3)«DELN(K>  

KK-K+3+NUMNP  • 

XDEL!I,IK+9)»DDEL!KK)  ' ' ; • 

DEL  T4( 1 , IK  + 9 l-DELN(KK) 

145 -CONTINUE  

150  CONTINUE  . 

C FIND  DISPLACEMENT#  8Y  AVERAGING*  AT  CENTER  OF  EACH  ELEMENT 

XDEL  !1,13)»!X0EL<  1#1)+XDEU2*1  )+XDEl<  3*  1)+ XDEL! 4*1)  )/4.00 

..: x DE  L ( 1*1 4 ) =■(  X D E L-M  » 2 > +X  D ELI  2*  2 ) + XDE  L-I  W-)  + XDE-L-1 4*-2-)4/  4-40 

XDEL l 1*15  ) • (XDEL ( 1,3) +XDEL (2*3) +XDEL! 3*3) +XDEL!  4*3) 1/4.00 
XDE L < 1* 16 )»!  XO EL ( 1*7) +XDEL (2*7 >+XDEL<  3*7) +XOEL  (4*7) 1/4.00 
XDEL  ( 1,  17  ) * ( XDEL(1,8)+X0EL!2»8)+XDEL13#8)+XD£L!4*8)  )/4.00 

- XDE  L I 1*  184  ■ ( X DE  L(-l*-9-M-X0  EL-1 2»-9.)-*XT)E  L43»  9 ) + X D £ L-(-4  *-944  A4-„0a 

D£LTA(1»13)-(0ELTA(1,1)+DELTA(2,1)+DELTA< 3* 1 > +DEL TA (4, 1 ) ) 
1/4.00 

DELTA!  1,14)  - ( DE LT A! 1* 2 > +DE LT A ! 2*  2)  + DF.LTA!  3*  2 )+DELTA  (4*  2 ) ) 

0 E L T aTi»‘1  5f* VdTl t‘aT r»T)  + DETta1  2 * 3 )+’dFlTM3 7F)+0 E L T A ( 4*7)7 

1/4.00 

i DELTA! 1, 16) » ! DELT A ! 1* 7 ) + DE LTA 1 2* 7 ) +DELTA!  3* 7 l+DELT A <4* 7) ) 

DELTA  Cl » 17)' •"(  DEL*  TaT  i'"»  sT*  OeXt  A (2*8  HOTlfATVisi +OEL  TA  (4*81) 
1/4,00 

DELTA! 1*16)«I DELTA !1» 9 ) +DE LTA I 2* 9 ) + DE LT A I 3*9 ) +DE LT A !4, 9 ) ) 

. - 1/4.00 

C CALCULATE  AVERAGE  STRAIN  •' 

DO  175  1-2,4  . . 

DO  175  J-13,18 

XDE.L  .{■!-*  J )j»XDEL  1-1,-i-l— 

DELTA! I, J)»DELTA!1,J) 

175  CONTINUE 

179  BSUM-0.00  ' ’ ; 

D 0- 1 8 0.  I - L,  6 ^ : : 

TSUM! I >-0.00 
180  CONTINUE 
CALL  OUAD 

DO-  .19.0.  I - 1,  6 

190  DEPS! I )=TSUM (I 1/0SUM 
C CALCULATE  STRESS  INCREMENT  . : 

DO  195  1-1,6  .-  * 

D S I G i I 4 » 0 .0.0 ; : 

00  195  J - 1 * 6 

195  DSIGID-  DSIG(I)  + CRZ!I,  J)*DEPS(J) 

CALL  YIELD. 

. C4LL  .MODI  F.Y 

200  CONTINUE 
RETURN 
END 
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I 


IROUTINE  TRISTF 


7A/7A 


OPT  = 1 


FTN  A.6+-A28 


7i 

“i 


i 


50 


SUBROUTINE  TR I S TF ( I I, J J ) 

INTEGER  CODE 

C0MMQN/8ASIC/VQL-,  NUMNP*  NUMEL»NUML A - 

COMMON/NATP/RO  ( 12  ) » E ( 9 » 12)  , E E ( 9 l 
COMMON /ARG/ X XX ( 101,  YYYUO  ).S(2A)*XX(  3), YY(3>, 

1CRZ(6,6>,  XI  (10)»S  IG(12)*N*M 

•COMMON  /N  P OA-T  A / X ( 200 ) * Y ( 200 ) * CODE ( 200 ) » NPNUM ( 10  »-20  )■■•-•• 

COMMON/ELDATA/3ETA(12>,  AL PHA ( 12 ) , TH ( 12 > , IX  ( 200, A ) , MATR IL ( 12 > 
COMM0N/TD/IMIN(100)*IMAX( 100) , JMIN( 25 ) * JMAX< 25)» MAX  I* MAX J*NMTL,NBC 
C0MM0N/8ASIC2/8ET»H*TIME*NLAY*IFLAG»IT»NTIME*LINC 

-COMMON / $ IG^-B£l)H*^-Sim  (6->-*-5i  G*A-(  12 *-50*6 It D$  IG+fa  )-*-H-6  GO  I 

COMMON/O I SP/0ELN1 (600 )* DEIN ( 600 >* DEL ( 600), GNM1 ( 600 > *GNM2( 600) 
COMMON /MASS/ A (600 ), 8(600) ,CM(8> 

COMMON/DELT/XDEU  A, 18) 

..-C  OM MON  / 0 E-L-T  R4-/-DS-G-T-M  A,  16  ) 

DIMENSION  0( 6,18) »C(18» 18) *XY( 9,9), SAVE (18,3) » BT ( 1 8 *6 ) * BTS I G (18)# 
1C3AR(10*18) *Z(10,10>*  CBD(10)»CC(10*A) *BLT (18*4)* 

28 AOO ( 1 8*  6 ) 

— D-I M E N S-I-0  N~&l-T-S4G-(  1 8->-*-&$-I  G44-8-) 

DIMENSION  DUM ( 10* 10 ) 

NLA«M 
NE  L*N 

— DO— 50 — I-=-l*  6 : = — 

SIG(I >=SIGMA(M*N* I) 


THI CK»TH ( NLA Y ) 
ISU3-  9 

XX-(-l-)=~X-XX-(  144- 
XX ( 2 ) “ XXX(JJ) 


xxx(isua) 

YYY(II) 
Y-Y-Y-(-J-J4- 


1*18) 


AO 

100 

150 


800 


X X ( 3 ) 

Y Y ( 1 ) * 

Y-Y-(-2-)-“- 

YY(3) * YYY< ISUB ) 

IF ( NTIME  * EQ .1 ) GO  TO  AO 
READ ( 2 ) ( (BT(I, J>* J=l*fa)*I 

GO— -TO- 801 

CONTINUE 
CALL  INTER 
CONTINUE 

00  —3  00— 14=1-*  6 

DO  300  J A =1 *1 8 
Q(  £ A*  J A)  = 0*  00 
0(1,2)=  XI(1)+  THICK 
-0-(-l,-U4»_X-U44-*-(.(-XHLCK)  **2)72^00-- 
3(3*13)=  X I ( 3)+  THICK 
0 ( 5 * 1 A ) = X I ( 2 ) + THICK 
Q ( 2 * fa ) ° 0(1*2) 

-Xt-3.,  1-64=— Q4-l-*-24— 

3(3*17)=  Q ( 5* 1A ) 

0 ( A*  3 ) = Q (1*  2) 

0 ( A * 5 ) * 0(1*2) 

3 ( 5. *.9  ) =_. Q_(  1*.2J 

0(5*13)=  0(1*2) 

0 ( fa  * 8 ) = 0(1*2) 

0(6*10)=  0(1*2) 

Q.(  2*l  5)  =-.  0Xl*-ll.).- 
0(4*12)=  0(1,11) 

0 ( A * 1 A ) = 0(1,11) 
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ROUTINE  TRISTF 


7*/7*  OPT  »1 


FTK  A.6+A28 


0 ( 5, 1 8 ) 3 0(1, 11) 
0(6,17)-  0(1,11) 
0 ( 5 * 1 5 ) 3 0(3, 18-1- 
•3(5,12)*  3(3,18) 
0(6,11)“  0(5,1*) 


THE  FOLLOWING  EXPRESSION  IS  ACTUALL  Y TM ICE  • THE  AREA-0F  THE  -T-RIAN6-LE-, 

AREA*  XX (l 
1XX(1')*YY(3 


AREA.3  XX  ( l J+YY(2)-XX(2)*YY(1)+XX(2)*YY(3)-XX(3)*YY(2)+XX(3)*YY(1)’ 


ZERO  C MATRIX 

DO  *00  I»l,18 

...  00  . *00- J«-l,18- 

*00  C(I,J)«  0.00 
FAC-l/AREA 
M 3 1 

N 3 1 

1 3 II 
J3  JJ 
GO  TO  250 
-.-7-5.  CON  TIN  UE- 

FAC»l/( AREA+THICK) 
M-10 
N3  1 0 

,.  I ■ 11  + *.. 

J3JJ+* 

FILL  C INVERSE  BY  PARTS 


250  C(M,N>"  (XXX(J)*  YYY ( ISUB )-  XXX(ISUB)*  YYY(J))*FAC 
C ( !“1 » N+  3 ) 3 (XXX(ISUB)*  YYY(  I ) - XXX  (I)*  YYY(  ISUB)  ) *FAC 
C ( M , N + 6 ) 3 ( XXX  ( I ) * YYY  ( J )-  XXX(J)+  YYY(I))*FAC 

. - C ( M +.1 , N.)  *.(  Y-Y  Y.l  J.).^-Y.Y-Y.l  I S UB-)  .MufcAC- : 

C(M+1,N+3)3(YYY(ISUB>-  YYY ( I ) ) *F AC  • 

C< M+l  ,N  + 8 ) 3 ( YYY < l )~  YYY(J))*FAC 
C(M  + 2,N)»(XXX(ISUB)-  XXX(J))*.FAC 

C ( A + 2,  N + 3.) 3 (XXX  UO^-X  XX  C IS  US ) ) *-F-AC: ™, — 

C ( M +2 » N + 6 ) 3 ( XXX  l J ) - XXX(I))+FAC 
C ( M+3, N+ 1 ) 3 C(M,N) 

C ( M + 3 , N+  * ) 3 C (M,N  + 3) 

,.C(  1+.l,.N.+  7.)3  .C.(M,M.+A) :_1_ 

C(M  + *,N  + 1 )=  C ( M+ 1 » N ) 

C (M  + *»  N + * ) 3 C (M+l,N+3) 

C(M+*,N+7)3  C (M+li N+6) 

C ( M + 5,.N+.l ) 3 ..  C ( M+.2  , N)._. 

C M + S,N  + * )-  C( K+2,N+3) 

C(M+5*N+7)=  C (M+2,N+6) 

C(M+6,N+2  )3-  C(M,N) 

C ( M.+  6,N.+.5.)  3 C.(  H , N.+.3.) :• 

C( M + 6, N + 8 ) 3 C (M,N+6 ) 

C(N  + 7,N+2)3  C ( M + 1 » N ) 

C ( M+  7, N*  5 ) - C ( M + l, N+3 ) 

.Cl  H+.7.,j!l  +.8.  ).= . ,.C..( M.+.1, Nib ) — _ 

C ( M + 8 , N + 2 ) 3 C ( M+2 , N ) » 

C(N+8,N+5)3  C ( M+2 , N+3 ) 
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^ROUTINE 


TRISTF 


74/74 


OPT  *1 


FTN  4*6+428 


C M + BfN+8)  ■ C ( M + 2 » N + 6 ) 

IF  ( ISU3.E0.10)  GO  TO  200 

I 503*-  '10  - -- 

GD  T3  75 
200  M*  Nt  A 
N * N E L 

FILL  XY  MATRIX 


500 


C 

C 

c 


00 
DO- 
X Y ( 
X Y ( 
X Y ( 

XY-4 

X Y ( 
X Y ( 
X Y { 

X Y ( 

X Y ( 
X Y ( 
X Y ( 

XY  ( 

XY( 
X Y C 
XY< 

— XY ( 

X Y ( 
X Y ( 
X Y ( 

— X Y ( 

X Y ( 
XY  ( 
X Y ( 

XX  t 

X Y { 
X Y ( 
X Y ( 
XY  ( 


600  I 
600 -J 
I J # J I 


1) 
2>« 
*3*)^ 
2 ) * 
3)* 
2>» 

3) * 

4 )  « 
7)* 
1)« 
*44* 

7 )  » 

n* 

4 ) * 
74  * 
5>* 

8 ) * 
, 6 ) ■ 
3^94* 
5# 5>* 

9 ) * 
6 J * 
-94* 

5) » 
8 ) c 
6 ) ■ 
94-« 


J*l#  9 

I-*l  *-9 

)*  0*00  . 

1*  00 
XXX  (I > 

-YY-Y4T4 

XXX  I J ) 

YYY  (J) 
XXXUO) 

-YYY  44-04 

XY  C X#  X ) 

XY ( 1 j 1 ) 

XY ( 1 # 1 ) 
-XY4  1#44— 
X Y ( 1 > 1 ) 
XY(hl) 

X Y ( 1#  1 ) 

X*44#44 

XY (1,2) 
XYU,2) 
XY<1#3) 

--X-Y  (4i3) 

XY  1 4#  2 ) 
XY<4,2) 

XY( 4, 3) 

X-Y-U*-34 

X Y ( 7#  2 ) 
XY(7,2) 

X Y ( 7#  3 ) 
„X*4-7,34- — 


COMPUTING  THE  LOWER  LEFT  QUADRANT  OF  THE  C MATRIX 


00-62  Q..XiuU,9 - 

DO  620  J*l#9 
OUM ( I, J ) *0* 00 
DO  620  K * 1 # 9 ■ 

oumxi.#  jx*  oumu-sx^vax  jju-tujujj- 

620  CONTINUE 

00  630  I » 10#  18 
DO  630  J n 1 # 9 

C U#J  >*0*00 - 

DO  630  < ° 1#  9 
9 n K +.  9 

C(I,J >«:(I»J)-C(i»K9)*0UM(K#J ) 

'530  ....CONTINUE- - 

: REARRANGE  THE  DISPLACEMENT  MATRIX 
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ROUTINE  TUSTF 


7A/74  OPT*  1 


FTN  <i.6+A20 


T 


v 

f - 


DO  730  I A*  1,  18 

DO  730  JA*.1,3 

J8*  JA+  6 

730  S AVE ( I A , J A ) » C ( I A , J B ) 

DO  740  IA*1,18 

• - . 00  740  J A*  7/  12  • 

J B*  J At  3 

740  C(IA,J4)»  C(IA,J9) 

DO  750  1 A*l, 18 

00-7-50 -44 j-1-,4 

JB-JAtl2 

750  C ( I A, JB) * SAVE ( 1 A > J A > 

C-ALCO  L ATE— THE—&Q— TR  AN-SPOSE-D— MATRIX- 

DO  300  J A* 1# 18 
00  800  I A*l, 6 

• - — BT  ( JA  » l A ). *4.  00 ^ 

00  800  MA* 1 * 18 

800  BT(JA,IA  > «*  BT  C JA, I A ) +C (M A, J A ) + Q ( IA,NA ) 
WRITE (2 ) ( (3TCI ,J ) >J*1.6)> 1*1, 18) 

-.801--.-  CONT  INUE-— 


CALCULATE  THE  BOT  * STRESS 
- DO  810— I- A*  1 >4 8 .... 


8TS IGC IA)*0. 00 
00  819  MA«1|6 

810  8TS IS  1 1 A ) “3TSI G(  I A ) + BT(  I A,  NA)  * SIG  (NA  ) 

ado  the" nonl'ine^Terhs  ~ : “ 


CALCULATE.  THE-C-BAR  .MATRIX-— 


IF  (NTIME.EQ.l)  GO  TO  815 

IF  ( (NTIME/20  J+20.EQ.  NTIME)  GO  TO  815 

REACT  34 — (-l.BLTC-I.j.J  L* 4®  1 , A 1,1*1,18.) 

GO  TO  661 
815  CONTINUE 

DO  020  I A * 1, 10 

03  82 0.-4-A*  1,18 - - . 

IB*  I A+7 

I F ( IA.EO .9)18*17 
I F ( IA.E3. 10)18*18 
320  ..C  B AR C I A^J  A.)  - ,C-(  IB, 4A.) .. 

CALCULATE  THE  Z BAR  MATRIX 

— 00  830  .14*1,10 

00  830  J A * 1 , 10 
830  Z ( I A, J A ) * 0.00 

Z(l,l)»  XI ( 1)*THICK 

Z 1 2.,2JL» Z.MjJ..) 

Z ( 3, 3 ) * Z(l,l) 

Z ( 6, 6 ) * Z ( 1 , 1 ) 
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iROUTINE  T31STF 


7A/7A  OPT-1 


FTN  A.6+A28 


-C- 

C 

C 


Z<3,A)»  XI  (2>*THICK 
Z ( A , 3 ) ® Z(3,A) 

Z<6,74«  Z (3,  A) 

Z ( 7»  6 ) ■ Z(3,A) 

U 9, 9)*  (THICK**3)/3+XI <1) 
Z(10,10>«  Z(9,9) 

Z ( 1,9 )» (THICK**2>72*XH1) 
2 ( 2> 1 0 » « Z(I*9) 

Z<9»1)«  Z(l*9) 

Z ( 10,  2 ) ° Z ( 1 , 9 ) 

Z ( 3,  5 ) « X-I-<-34*-T4«G-K 

ZC5,3)»  Z(3,5) 

Z<6,8)»  Z ( 3, 5 ) 

Z ( 8, 6 ) ° Z ( 3, 5 ) 

-Z-(  4 , 5-)  *~T+»-I-CK- 

Z ( 7, 7 ) » Z (<,><* ) 

Z(5,5)»  X I (6  )*  THICK 
Z ( 8 , 6 ) * Z(5,5) 

- Z-(  A,  5 ) »~-X-tm-*-THKK^ 
Z<5,A)»  Z(A,5) 

Z ( 7 » 8 ) » Z<A,5> 

Z ( 8, 7 ) “ Z<A,5) 


CALCULATING  THE  CC  MATRIX 

DO  8A0  I A“1 , 10 
.....;..C&D-MA4-»-0.-&&- 


00  8A0  J A" 1, 18 
8 AO  CBD(IA)-  CBD  ( IA  )+•  CB'AR  ( I'A,  J A )*  DELTAtII,JA) 
00  850  IA=1,10 

0O~85O~M«4*A — 

850  CCCIA, JA)»0.00 
CC.(1,1)«  CBO  ( 1 ) 

CC(2»A)»  CBD ( 1 ) 

C C ( l-f-AT*--  CB-D42-) 


855 


CC ( 2*  2 ) ® CBO ( 2 ) 

CC ( 3, 3>»  CBD (3  ) 

CC ( A, 3 ) ■ CBD (A ) 

- C C-t-5,  3 )-■*-•  -C&W--54- 
CC(6,3)»  C 00  16  > 

CC ( 7, 3 ) ■ CBD (7 ) 

CC 19, 3)«  CBD ( 8) 

-C  C (9,4r)-»~C-B-D49  >— 

CC ( 10, A ) ® CB D ( 9 ) 

C C ( 9,  A ) ® CBDMO) 

CC ( 10,2)®  CBD(10> 

— 0-0  -8-5  5 — I * 1-f  l-O  - 
DO  855  J«1,A 
DUMt I, J J-0.00 
DO  855  <®1,10 

0U1.(  I.,  J-).»OUM(.I.,.J-) -+Z4-I,K.).*CC  CK,  J ) 
CONTINUE 
DO  860  I <*  1, 18 
DO  ?60  J«1,A 

3 L T ( I ^J-).«.0-.00- 

00  660  K-  1, 10 
Q.LT(I,J)*BLT(I,J)+CBAR(KjI)+DUM(Kj»J) 
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ROUTINE  TRISTF 


7W7A  OPT “1 


FTN  A. 6 + A28 


860  CONTINUE 

WRITE  (3)  UBLTCIf  J)fJ**l,  A),  1-1,18) 

061  CONTINUE  - 

00  070  I A - 1 > 10 
BLTSIGUA)*  0,00 
DO  870  JA»ljA 

370  BLTSIG(IA)-  BL TS I G< I A ) ♦ Bl T ( I A, J A >♦  S IG ( J A ) 

00  830  I A»  1>  18 
BSIG(IA)-  BrSIG(IAJ+  BLTSIG(IA) 

380  CONTINUE 

: AOD I NG  -THE  --T-R-I A N G U L A R - B-5  1-6-HAT R I-GE-S-T-0  F-ORM-  T HE-OU AOR^A-T-ERAV-MATR^X- 

I F (I FL  AG • NE • 2 ) GO  TO  1600 
: EVALUATE  STRAIN  INCREMENT 

BSUM-9$UM+AREA/2.0*THICK 

DO-  - 5 2 5>  - -I  o-l  >-1-8 « t : 

00  525  J*1,A 

BAOOII, J )-  BT(IiJ)i  BLT(I,  J > 

525  CONTINUE 

DO  526-1  »1#  1-8 : 

DO  526  J-  5, 6 

526  BAOOdiJ  )»  BTd,J) 

00  550  I » 1 1 6 

- 00-5  50-J*-l-f-ia- 


550  TSUMm*TSUMtI)+BADDt  Jjl  )*XDEl  (II* J ) 

GO  TO  1000 
1600  CONTINUE 

— K«..«3+II»2 

L«  3+JJ-2 
DO  910  MA  * If  2 
MJ«  <MA-1>*  6 

.N  J tt  A.* U.A — 12- 

DO  920  I A — 1 # 3 
JA«  lArl 

$(K  + NJ*  J A ) ■ S ( K+  NJ+  J A ) + BSIGCMJ+*IA> 

.920  S ( L + N J-+  .J.Al.p--S.a-+-N.J>-JA >..+  -BSIG  <4W-+_tA*^34 - 

DO  910  I A»  If  3 
DO  910  J A * 1 f A 
IJ-(JA-l)+3 

910  »S  ( N J + I J*- 1 A ).*— S-(  J-*l  A-4-+-  B S-I-G  C-12 .+-(  MA^Jr-)  niI  A4  /4*4 

I F ( NTIM6.GT.1 ) GO  TO  1000 
C ASSEMBLE  ELEMENT  MASS*  MATRIX 

DEN«RO(MATRIL(NLAY) )+ THICK 

XMA  S S.“  AR,E.A£D.EN  /.8  *00- _ - 

CMC  II  ) * C <M  ( II  ) + XMASS 
CMC  JJ) -CM  < JJ  HXMASS 
CMC  II+A)»CM< II ) 

CM!  J J .f  A.UC  il  l JJ.1- „ 

1000  CONTINUE 
RETURN 
END 


IB  OL  I C ‘R  EFERENCE  map  Tr  ■ 1 ) 
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BROUTINE  YIELD 


74/74 


4*6+428 


SUBROUTINE  YIELO 
INTEGER  CODE 

COMMON /BASIC/ VOL*  NUMNR*  NUM€l*  NUMLA  

C0MM3N/MATP/R0(12)*E(9,12  *§E(? 

COMMON/ ARG/ XXX  ( 10)*YYY(10)»S(24)»XX(3)»YY(3)* 

^COMMON /NP DATA  / X(  200)  »Y  I 200  ) *COOE  ( 200) » NPNUM( 10*  20) — — — 

CaMM0N/ELDATA/BETA(12)*ALPHA< 12)»TH(12>* IX( 200*  A )*MATRIL(12) 

C0MM3N/RESULT/0(6*6)*C(6*6j^<^S(6*6) 

COMMON /TD/ 1 MINI  100)* IMA  X ( 100 ) * JM I N(25 ) » JM AX ( 25) * MAX  I* MAX J » NMTl  i 

.....  COMMON/BA  SIC  2 AB  E T > H »-T I M E *-N4.-A-Y-*.I  FLAG*-  I-J-+-NT  I M E*4-IN  C ■ 

COMMON / S I GM / 8 SUM* TSUM ( 6 ) * S I GMA ( 12  * 50* 6 ) * D S IG ( 6) * F ( 600 ) 

COMMON/DI SP / DE  LN1 (600 ) » D ELN ( 600 ) * DEL ( 600) » GNM1 (600 ) #GNM2( 600 ) 
COMM0N/MASS/AI6O0 )>  B(600)*CM(8) 

COMMON / P L Y LD- AS-I-G-Y-MB-)  > D EP  S U>-> — — 

DIMENSION  SIGT( 6) * SIGND ( 6 ) * S I GN ID ( 6 ) » S I GO ( 6) 

DO  15  1=1*6 

SIGT(I)aSIGMA(M*N*I)+DSIG(I > 

'calculate'  theTi i Vi atorTc“c o m> o nentTof' ‘ st r es $ 

$KNNlaSIGT(l)+SIGT(2>+SIGT(3) 

DO  30  1=1*3 

SI  G D.(  I)  =S-UH4-I-USKNN-L/  3-*  0 

CONTINUE 
DO  31  Ia  4*6 
S I3D( I ) «SIGT  (I > 

...  CONTINU  E ; 

YLDT»FSIGD(l)°*2  + SIGD(2)+<'2  +SIGD(3)**2  + 2.  *(  SIGD  (4)  **Z  ♦ 

1 SIGD(5  )*  + 2 + S IGD ( 6 ) + + 2 ) -(2./3.)  *SIGY (M  ) + * 2 

n'AS  T I C ‘'FLOW  aHALSG  OC'C  UR  RE  oT(T'L  A Y E.R  = M*Tl¥mFnT  = N 

CORRECT  STRESSES  FOR  PLASTICITY 
CALCULATE  LAMBDA-BAR  ^ ^ 

AO « SI  GO ( ll + * 2 + S I GO T 2 > * * 2 + SI G '67  3 ) *+Z  + 2 . * (' sTgTTa  ) * + 2 

lB0-  SIGT?i?+SIG0(l?+$IGT*2)*$IGD(2)+3IGT(3)*|IGD(3)+2»*( 

1 SIGH-4 ) *.S  I£a(Al+S  I6X  (-5  ■)  ASJ  G04-54  +-S4£XiAl-»-S  1 GO  ( 6-)-)— 

CO®  Y L D 

DISC  = 3 Q + * 2 -AO+CQ 
IFIOISC.LT.  0.0)  OISC=0.00 

,00  XLMSTR-  £0./  (B0.+  .SQR.UD.LSC)  ) — 

CALCULATE  PLASTIC  STRESS  INCREMENT 
DO  AO  1=1*6 

DS I G( 1 ) =D SI G ( I )-S ISD( I) +XL MS TR 

) CONTINUE — ; 

)C  CONTINUE 
)0  CONTINUE 

RETURN  


+2  ♦*( SIGD (A ) * + 2 
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